Home  Contact Us  Orders: US  Orders: other countries   Review Shopping Cart

Serious mathematics, written with the reader in mind.
[MsE] Matrix Editions

 

 

Monte Carlo Program

This program requires a Pascal compiler.

The nine lines beginning function Rand and ending end are a random number generator.

The six lines beginning function randomfunction  and ending end define a random function that gives the absolute value of the determinant of a 2 x 2 matric with entries a,b,c,d.  For an n x n matrix, you would enter n squared "seeds".  You can name them what you like; if n=3, you could call them x1, x2, ..., x9 instead of a, b,...i.  In that case you would write x1:=rand(seed); x2:=rand(seed), and so on.  To define the random function you would use the forrmula for the determinant of a 3 x 2 matrix.

Here is the program. After the program we give an example showing how to use it.




program  montecarlo;

const lengthofrun = 100000;


 var S,V,x,intguess,varguess,stddev,squarerootlerun,
 errfifty,errninety,errninetyfive:longreal;
 i,seed,answer:longint;

function Rand (var Seed: longint): real;
{Generate pseudo random number between 0 and 1}
const  Modulus = 65536; 
     Multiplier = 25173;
     Increment = 13849;
begin
 Seed := ((Multiplier * Seed) + Increment)  mod Modulus;
 Rand := Seed / Modulus
end;

function randomfunction:real;
var a,b,c,d:real;
begin
 a:=rand(seed);b:=rand(seed);c:=rand(seed);d:=rand(seed);
 randomfunction:=abs(a*d-b*c);
end;


begin
Seed := 21877;
repeat
S:=0;V:=0;
for i:=1 to lengthofrun do 
begin
      x:=randomfunction;
      S:=S+x;V:=V+sqr(x);
 end;        
intguess:=S/lengthofrun; varguess:=V/lengthofrun-sqr(intguess);
stddev:= sqrt(varguess); 
squarerootlerun:=sqrt(lengthofrun); 
errfifty:= 0.6745*stddev/squarerootlerun;
errninety:= 1.645*stddev/squarerootlerun;
errninetyfive:= 1.960*stddev/squarerootlerun;
writeln('average for this run = ',intguess);
writeln('estimated standard deviation = ',stddev);
writeln('with probability 50%, your error is less than
',errfifty);
 writeln('with probability 90%, your error is less than
',errninety);
 writeln('with probability 95%, your error is less than
',errninetyfive);
writeln('another run? 1 with new seed, 2 without');
readln(answer);
 if (answer=1) then
 begin
 writeln('enter a new seed, which should be an integer');
 readln(seed);
 end;
until (answer=0);

end.

Example of using the Monte Carlo program:

Note: In Pascal, x^2 is written sqr(x).

To compute the area inside the unit square and above the parabola y=x^2,  you would type

function functiomfunction:real;
var x,y:real;
begin
x:=rand(seed);y:=rand(seed);
if (y-sqr(x) <0 ) then randomfunction:=0
else randomfunction:=1;
end

 

Back to Matrix Editions home

 

 

 

 

 

Last updated March 12, 2003