Approximating an integral using Monte Carlo Method

Question

I wrote a solution for

Calculate the value of the integral I = $\int_0^\pi sin^2(x)dx$ using the Monte Carlo Method (by generating $ 10^4 $ uniform random numbers within domain [0, π] × [0, 1]). Do the calculation using LCG (Linear Congruent Generator) with a = $7^5$ , c = 0 and m = $2^{31}$ − 1.

In my solution I use the formula for LCG, but my random values keep getting smaller and smaller each time. Also, no initial value is mentioned and if I used 0, then the random numbers would always be 0, so I used r=1.

r = 1
hit = 0
for i = 1:1000
  r = (mod(7^5*r + 0, 2^31-1)/(2^31-1))
  x = r * pi
  r = (mod(7^5*r + 0, 2^31-1)/(2^31-1))
  y = r
  if (y <= sin(x)*sin(x))
    hit = hit + 1
  end
end
hit/10000 * pi

I get 3.1322 as a result, but the actual value of the integral is approx. 1.57. I don't see how can I get anything bigger than 1 if I divide hit/total, though that's what my notes suggest me to do...

Answer 1

Division by $2^{31}-1$ messes up the LCG. Keep r integer; division is needed to get x and y from it.

r = mod(7^5*r + 0, 2^31-1)
x = r*pi/(2^31-1)
r = mod(7^5*r + 0, 2^31-1)
y = r/(2^31-1)

Answer 2

The formula for the LCG is $r_n= mod(a \, r_{n-1} + c, m)$. It's ok that you then normalize the result it , diving it by $m$, but you should not store the result in $r_n$, but in a separate variable. The sequence of numbers that the LCG generates are integers.