Pick random numbers between 0 and 1, evaluate the function and average
the result. This is simple
The difference between any function of N statistically independent points and its true value is due to the errors introduced by the points.
The difference squared is
In addition, its expectation value is
The uncorrelated values for i ¹ j are always assumed to be zero in the statistical ensemble. It is this assumption that enables us to arrive at a reasonable approximation to the ensemble average using a single element of the ensemble.
The term in the sum in can be found by leaving out the i'th data point. That is
Equation is a sum over the result of omitting each term in turn. If a minimization method such as the Amoeba .htm or Nlfit .htm was used to find f it needs to be rerun with the point omitted. If Lagrange .htm interpolation was used to find f, the interpolation needs to be made again with the point omitted. Since the loss of a single point does not move Nlfit or the Amoeba far from equilibrium, the new value of f – with a missing point - can usually be found with only a few steps from the original – with all points. In the case of the Amoeba the full output array should be used as input for each to the new sets. Convergence should be relatively rapid since one is nearly at the minimum in each case. As long as special points such as the ends are removed, it is always possible to "randomly" select M of the points and then to multiply by N/M to arrive at the final error estimate.
Psuedo Random numbers
Random numbers are not wanted. What is desired is a set of completely reproducible numbers that are the same on all computers and in both C and Fortran. In addition, it must be true that the numbers give all averages the same as would be given by actual random numbers. In practice the numbers should be re-initialized with the numbers needed to recreate the desired set written out as
17894561 3678991
18345899 234567899
...
In this way, when the code bombs after 23 hours of running, the situation just before the bomb can be recreated by simply entering these numbers.
With this many random numbers, the "randomness" of the set can be a problem for some calculations[1]. The method used here, first suggested by MacLaren and Marsaglia[2], is to insert numbers into a table with a multiplicative congruential method with one multiplier and seed and to extract them from the table with a multiplicative congruential method with a different multiplier and seed. random.for contains the common
COMMON/RAN/IX,IY,ITAB(128)
The file random.h contains
struct {int ix,iy,itab[128];} ran;
int rseed(int ix,int iy);
double rndmf();
The seeds are ix and iy. The random number set is initialized by calling rseed(ix,iy). This call causes 128 numbers to be fed into the table itab[128] with one random number generator. The call x=rndmf() gives a random number base on ix,iy and itab, between 0 and 1. It selects this by advancing ix and iy, pulling a value from a randomly selected table entry and replacing that with a new value.
This can be re-initialized by printing out ix and iy and then calling rseed. This wastes 128 random numbers, but gives a completely determined new set of random numbers.
Estimate
Do this both by choosing random numbers and by using a mid-point grid in all six dimensions. Make w small enough that the integral is w6p3. Send me the code, the answers with errors for the Monte-Carlo results and some discussion. ../integration/MonteCarlo/FMonte/Fmonte.doc#Biasedasrandom .htm
Rgauss=-6;
For(i=0;i<12;++i){Rgauss+=rndmf();}
Thus
or
This means that the above sum is 0 +- 1. If I want to produce a peak of height P = 36 +- 6 with a Gaussian distribution on the six. I simply take
[1] R. L. Coldwell, "Correlational Defects in the Standard IBM 360 Random Number Generator and the Classical Ideal Gas Correlation Function", J. Computational Phys. 14, 223 (1974).
[2] M. D. MacLaren and G. Marsaglia, J. Assoc. Comput. Mach. 12 (1965),83
[3]
J. M Hammersley and D. C. Handscomb, Monte Carlo Methods, John Wiley (1964) p.39