Here is a console mode UVS program which generates some normally distributed numbers, plots them in a histogram, and tests the mean and variance to see if they are what they are supposed to be:
// uvs_example.cpp: example usage of uvs library #include <iostream> // for std::cout #include "uvs/uvs.h" // for uvs library
We include iostream because we are going to refer to std::cout later. The other include brings in the entire uvs library.
void
Monte_Carlo()
{
// process generator
typedef uvs::Good1 RNG;
RNG rng; // the random number generator
double mu1 = 1.25; // the population mean
double sig1 = 0.75; // the population standard deviation
uvs::NormalPG<RNG> pg1(&rng, mu1, sig1);
The function "Monte_Carlo" will do the things alluded to above. Here we have chosen to use Good1 for our process generator. We define it indirectly using a typedef to make it easy to change this decision later if we want to. We are going to generate normal numbers with a mean of 1.25 and a standard deviation of 0.75. Thus we declare a NormalPG<RNG> called pg1 which keeps a pointer to the random number generator and is initialized with the selected mean and standard deviation. This object is an example of a "smart process generator" -- one which knows many useful things about the distribution it generates. We will take advantage of some of these capabilities when we do our hypothesis testing below.
// sample statistic double samp1_min = -1.0; // histogram lower limit double samp1_max = 3.5; // histogram upper limit int samp1_nbin = 20; // histogram bin count uvs::CoHistogram samp1(samp1_min, samp1_max, samp1_nbin);
The sample statistic of choice is a histogram for continuous variables called CoHistogram. We initialize it with some bin parameters which give the lower and upper limit of the histogram and the number of internal bins. This class automatically adds one bin on each side to hold the values which fell below or above the limits. In this case we can put in the mean minus and plus three sigma since we know the distribution we are expecting. In more realistic examples it is sometimes necessary to experiment with these parameters. Two things you should know:
// generate and show the samples int npass = 10000; // number of samples to generate for (int jpass = 0; jpass < npass; ++jpass) samp1(pg1()); // generate and record a sample samp1.show(std::cout); // show the histogram as text
This uses pg1 (the normal process generator) to generate the random variables and send them straight to samp1 (the histogram). The call to samp1.show creates a text representation of the histogram, which looks like this:
sample size = 10000
mean = 1.24299
standard deviation = 0.750043
variance = 0.562509
smallest = -1.79845
largest = 4.05354
x % cum % 0% 5% 10% 15% 20% 25% 30% 35% 40% 45% 50%
----------- ----- ------ |....|....|....|....|....|....|....|....|....|....|
-1.000 0.17 0.17 :
-0.775 0.19 0.36 :
-0.550 0.45 0.81 :
-0.325 1.01 1.82 :*
-0.100 1.90 3.72 :**
0.125 3.11 6.83 :***
0.350 4.75 11.58 :*****
0.575 7.18 18.76 :*******
0.800 9.53 28.29 :**********
1.025 9.87 38.15 :**********
1.250 11.79 49.94 :************
1.475 12.03 61.97 :************
1.700 10.91 72.87 :***********
1.925 9.09 81.97 :*********
2.150 6.70 88.66 :*******
2.375 4.51 93.18 :*****
2.600 3.47 96.65 :***
2.825 1.75 98.40 :**
3.050 0.90 99.30 :*
3.275 0.44 99.74 :
3.500 0.16 99.90 :
0.10 100.00 :
Note that the mean was supposed to be 1.25 but was actually 1.24299. Pretty close, but could that really be statistical sampling error? That is what hypothesis testing can tell us. Here is what it looks like:
// test mean double alpha = 0.01; // probability of rejecting if null hypothesis is true uvs::HypMean hmean(samp1.xbar(), samp1.s(), samp1.n(), pg1.mean()); // the hypothesis bool reject_mean = uvs::reject_two_tail(hmean, alpha); // the test std::cout << "reject mean: " << reject_mean << '\n';
We define a risk we are willing to take that the following test will say the mean was wrong even when it was actually right. This probability is called "alpha". Because pg1 is a smart process generator we can ask it what the mean of its output should be. We will test the hypothesis that the actual mean is what the distribution expects. (We might also have used the parameter mu1, but this would only work for the normal distribution.) Then we test the hypothesis using a two tail test (because we don't care if the mean is too low or too high) using the alpha we decided on. We expect the value printed out to be 0 (i.e. false) and so it is.
We can do a similar test on the variance (which is equivalent to testing the standard deviation) as follows:
// test variance uvs::HypVariance hvar(samp1.s(), samp1.n(), pg1.stdev()); // the hypothesis bool reject_variance = uvs::reject_two_tail(hvar, alpha); // the test std::cout << "reject variance: " << reject_variance << '\n'; }
The hypothesis is different from the last example but the alpha value and the two tailed test are the same. As expected, the difference between the sample variance and what the process generator says it should be is not significant.
Now all we have to do is to call the function:
int
main()
{
Monte_Carlo();
return 0;
}
This example program is included in the uvs distribution.