Simple Random Number Generation

using System;

namespace TestSimpleRNG
{
    /// <summary>
    /// Test application for the SimpleRNG random number generator.
    /// This verifies that the random numbers have the expected 
    /// distribution using a standard statistical test.
    /// Unfortunately the test is more complicated than the generator itself.
    /// 
    /// For more information on testing random number generators, see
    /// chapter 10 of Beautiful Testing by Tim Riley and Adam Goucher.
    /// </summary>
    class Testharness
    {
        /// <summary>

        /// </summary>

        static void Main(string[] args)
        {
            // Test the core uniform generator using the Kolmogorov-Smirnov statistical test.
            KSTest();

            // Test the derivted distributions by looking at their means and variances.
            TestDistributions();
        }

        static void KSTest()
        {
            /// Kolmogorov-Smirnov test for distributions.  See Knuth volume 2, page 48-51 (third edition).
            /// This test should *fail* on average one time in 1000 runs.
            /// That's life with random number generators: if the test passed all the time, 
            /// the source wouldn't be random enough!  If the test were to fail more frequently,
            /// the most likely explanation would be a bug in the code.
            
            SimpleRNG.SetSeedFromSystemTime();

            int numReps = 1000;
            double failureProbability = 0.001; // probability of test failing with normal input
            int j;
            double[] samples = new double[numReps];

            for (j = 0; j != numReps; ++j)
                samples[j] = SimpleRNG.GetUniform();

            System.Array.Sort(samples);

            double CDF;
            double temp;
            int j_minus = 0, j_plus = 0;
            double K_plus = -double.MaxValue;
            double K_minus = -double.MaxValue;

            for (j = 0; j != numReps; ++j)
            {
                CDF = samples[j];
                temp = (j + 1.0) / numReps - CDF;
                if (K_plus < temp)
                {
                    K_plus = temp;
                    j_plus = j;
                }
                temp = CDF - (j + 0.0) / numReps;
                if (K_minus < temp)
                {
                    K_minus = temp;
                    j_minus = j;
                }
            }

            double sqrtNumReps = Math.Sqrt((double)numReps);
            K_plus *= sqrtNumReps;
            K_minus *= sqrtNumReps;

            // We divide the failure probability by four because we have four tests:
            // left and right tests for K+ and K-.
            double p_low = 0.25 * failureProbability;
            double p_high = 1.0 - 0.25 * failureProbability;
            double cutoff_low = Math.Sqrt(0.5 * Math.Log(1.0 / (1.0 - p_low))) - 1.0 / (6.0 * sqrtNumReps);
            double cutoff_high = Math.Sqrt(0.5 * Math.Log(1.0 / (1.0 - p_high))) - 1.0 / (6.0 * sqrtNumReps);

            Console.WriteLine("\n\nTesting the random number distribution");
            Console.WriteLine("using the Kolmogorov-Smirnov (KS) test.\n");

            Console.WriteLine("K+ statistic: {0}", K_plus);
            Console.WriteLine("K+ statistic: {0}", K_minus);
            Console.WriteLine("Acceptable interval: [{0}, {1}]", cutoff_low, cutoff_high);
            Console.WriteLine("K+ max at {0} {1}", j_plus, samples[j_plus]);
            Console.WriteLine("K- max at {0} {1}", j_minus, samples[j_minus]);

            if (cutoff_low <= K_plus && K_plus <= cutoff_high && cutoff_low <= K_minus && K_minus <= cutoff_high)
                Console.WriteLine("\nKS test passed\n");
            else
                Console.WriteLine("\nKS test failed\n");
        }

        // Convenience function for TestDistributions()
        static void PrintResults
        (
            string name, 
            double expectedMean, 
            double expectedVariance, 
            double computedMean, 
            double computedVariance
        )
        {
            Console.WriteLine("Testing {0}", name);
            Console.WriteLine("Expected mean:     {0}, computed mean:     {1}", expectedMean, computedMean);
            Console.WriteLine("Expected variance: {0}, computed variance: {1}", expectedVariance, computedVariance);
            Console.WriteLine("");
        }

        // Verify that distributions have the correct mean and variance.
        // Note that sample mean and sample variance will not exactly match the expected mean and variance.
        static void TestDistributions()
        {
            const int numSamples = 100000;
            double mean, variance, stdev, shape, scale, degreesOfFreedom;
            RunningStat rs = new RunningStat();

            // Gamma distribution
            rs.Clear();
            shape = 10; scale = 2;
            for (int i = 0; i < numSamples; ++i)
                rs.Push(SimpleRNG.GetGamma(shape, scale));
            PrintResults("gamma", shape*scale, shape*scale*scale, rs.Mean(), rs.Variance());

            // Normal distribution
            rs.Clear();
            mean = 2; stdev = 5;
            for (int i = 0; i < numSamples; ++i)
                rs.Push(SimpleRNG.GetNormal(2, 5));
            PrintResults("normal", mean, stdev*stdev, rs.Mean(), rs.Variance());

            // Student t distribution
            rs.Clear();
            degreesOfFreedom = 6;
            for (int i = 0; i < numSamples; ++i)
                rs.Push(SimpleRNG.GetStudentT(6));
            PrintResults("Student t", 0, degreesOfFreedom / (degreesOfFreedom - 2.0), rs.Mean(), rs.Variance());

            // Weibull distribution
            rs.Clear();
            shape = 2; scale = 3;
            mean = 3*Math.Sqrt(Math.PI)/2;
            variance = 9*(1 - Math.PI/4);
            for (int i = 0; i < numSamples; ++i)
                rs.Push(SimpleRNG.GetWeibull(shape, scale));
            PrintResults("Weibull", mean, variance, rs.Mean(), rs.Variance());

            // Beta distribution
            rs.Clear();
            double a = 7, b = 2;
            mean = a / (a + b);
            variance = mean*(1 - mean) / (a + b + 1);
            for (int i = 0; i < numSamples; ++i)
                rs.Push(SimpleRNG.GetBeta(a, b));
            PrintResults("Beta", mean, variance, rs.Mean(), rs.Variance());
        }
    }

    class RunningStat
    {
        int m_n;
        double m_oldM, m_newM, m_oldS, m_newS;
 
        public RunningStat()
        {
            m_n = 0;
        }

        public void Clear()
        {
            m_n = 0;
        }

        public void Push(double x)
        {
            m_n++;

            // See Knuth TAOCP vol 2, 3rd edition, page 232
            if (m_n == 1)
            {
                m_oldM = m_newM = x;
                m_oldS = 0.0;
            }
            else
            {
                m_newM = m_oldM + (x - m_oldM)/m_n;
                m_newS = m_oldS + (x - m_oldM)*(x - m_newM);
    
                // set up for next iteration
                m_oldM = m_newM; 
                m_oldS = m_newS;
            }
        }

        public int NumDataValues()
        {
            return m_n;
        }

        public double Mean()
        {
            return (m_n > 0) ? m_newM : 0.0;
        }

        public double Variance()
        {
            return ( (m_n > 1) ? m_newS/(m_n - 1) : 0.0 );
        }

        public double StandardDeviation()
        {
            return Math.Sqrt( Variance() );
        }
    }

}