November 4, 2014 --- Class 20 --- Testing Hit or Miss Method (last step)
Sample Mean Method, Inverse Transform Method
Testing Hit-or-Miss Method (last step)
In our last class we learned about the hit or miss method and ran
a code pi3 to test it. You wrote a script, hopefully similar to this:
[sw246-00:~/src/misc] sg% mo ~/runpi3.csh
#!/bin/csh
foreach num (10 50 100 250 1000 2500 5000 10000)
pi3 $num 10000 |aver |awk '{print $1,$2,100*$3}'
end
The error coming out of aver is multiplied by 100 because
we want the width of the distribution of pi3, output,
not the error in the mean of the 10000 trials in each run of pi3.
A typical run looks like this:
[sw246-01:~] sg% runpi3.csh
1.0000000000e+01 3.1430400000e+00 0.518381
5.0000000000e+01 3.1412960000e+00 0.230508
1.0000000000e+02 3.1401240000e+00 0.166837
2.5000000000e+02 3.1409952000e+00 0.103798
1.0000000000e+03 3.1424588000e+00 0.0522629
2.5000000000e+03 3.1413792000e+00 0.0330492
5.0000000000e+03 3.1419865600e+00 0.0229892
1.0000000000e+04 3.1417096800e+00 0.0162435
How Big Should the Errors Be?
In the last class we introduced a generating function and
almost completed the analysis to find the variance of the
number of hits. Here is the result:
F_N(p,q) = (p+q)^N
If we want to find the average number of hits = *, we can
apply the operator p X partial derivative w.r.t. p
We can also get ** by applying that operator twice. After applying
the differential operators, we can set p+q = 1
We find ** = N p and ** = N p + N(N-1) p^2 . So the variance
** - **^2 = N p (1-p)
Our pi3 program takes i and multiplies it by 4/N to get a number
that should be Pi. Taking these factors into account, the error
in pi, proportional to the square root of the variance should be
4 \sqrt( p(1-p)/N )
Using the fact that p=pi/4 above, we found that
the numerical result agrees prefectly with the prediction based on the
generating function F_N(p,q).
You can also write an awk script to print the absolute value of the
error for each trial. Here is the awk commaand:
{ err=3.1415926-$2; if(err<0) err= (-err); print $1,err}
Now you many apply this command to your file before sending it to
aver and then plot the error vs size. Use logarithmic scale for
the x and y axes. Can you verify that the error falls like
1/sqrt(size) ? We did not do this in class.
Sample Mean Method
The sample mean method is based on the theorem of calculus that relates
the integral to the mean value of the function. It is a lot like the
rectangle rule, but instead of taking uniform values in x, we take
random values in the range of x. The error only falls like 1/sqrt(N),
so it is a very slow method and not generally suitable for low
dimensional integrals. Section 11.2 on pages 421--424 discusses
the sample mean method and Problem 11.7 (pg 424) is similar how
we test this method in class.
Section 11.4 on pages 426--429 is important and I have been known to
assign problem 11.10, but I don't think I will do it this year. You
should consider doing it on your own. Let me know if you have questions.
We ran the sample mean method for the same case as the for the hit or
miss method. The class was given 10 minutes to produce a shell script
and run pi2, which implements sample mean. We found that the error is
smaller with sample mean. On the board I showed that the sample mean
error falls like 1/Sqrt[n] were n is the number of points sampled. The
prefactor is the Sqrt[ variance of the function over the range of integration].
I got a little messed up demonstrating that.
Here is a sample run with pi2 replacing pi3:
1.0000000000e+01 3.1397844922e+00 0.282752
2.0000000000e+01 3.1400366551e+00 0.199547
5.0000000000e+01 3.1389893764e+00 0.126696
1.0000000000e+02 3.1430579796e+00 0.0891922
2.0000000000e+02 3.1422493818e+00 0.0624787
5.0000000000e+02 3.1418025705e+00 0.0401641
1.0000000000e+03 3.1422376989e+00 0.028339
2.0000000000e+03 3.1415877008e+00 0.0199763
5.0000000000e+03 3.1415294028e+00 0.0125688
1.0000000000e+04 3.1414853823e+00 0.00897418
The output of pi2 is 4 * \int_0^1 dx \sqrt(1-x^2)
I forgot about the factor of 4 in class. The variance of the
integrand is 2/3-(pi/4)^2 = 0.0498164. We need to take the square root of
variance and then multiply it by 4 because pi2 multiplies the integral
by 4. The result is 0.892784, so the error in the sample mean method should
be 0.892/Sqrt[n]. Looking at the table above, for n=100, we have 0.08919,
which agrees very well with expectation. Life is good!
We talked about getting the error down by multiplying and dividing the
integrand by a probability density close in shape. If we then sample
according to the probabilty density, the variance of the integrand/density
may be smaller than the variance of the integrand itself. Thus, it is
useful to be able to create random numbers with non-uniform probability density.
Inverse Transform Method
Sometimes we are interested in generating random numbers with a non
uniform probability density. There are a number of ways of doing
this and we shall study a few. This is covered in CSM section
11.5, starting on page 429. The first method we will study is
called inverse transform. You may be used to the idea of the
probability density function. For a real variable, we have
$\int_{-\Infinity}^{+\Infinity} p(x) dx = 1 $.
We can define the monotonically increasing function called the
cumulative probability distribution by
$ P(x) = \int_{-\Infinity}^{x} p(x') dx' $.
Of course, $ d P(x)/dx = p(x)$. P(x) is a number between 0 and 1.
If we can solve the equation $ P(x)=r$ for $x(r)$, then if we
generate uniform random numbers r, then the numbers x(r) have the
distribution we want. So, the key is we start with p(x), and
the method only is useful if we can find P(x) and then solve
P(x) = r
for x as a function of r.
A very simple example is a constant probablility for x to be between
two numbers a and b. In this case, we easily showed that
P(x) = (x-a)/(b-a) for ab. Solving P(x)=r, we find
x = a + (b-a)r
This simple relation is exactly what one would intitutively select for
generating numbers uniformly distributed between a and b.
We will study other distributions in our next class.
We are planning a makeup class on Thursday, November 13 at 6 p.m.
*