Monte Carlo simulation using LCG

[kisslaszlo123]

Hello,

I try to code a VBA Monte Carlo simulation that calculates the value of Pi with the help of a linear congruential generator. However something is wrong with my sub. Can someone please take a look and explain me what is the problem here?

Option Explicit

Sub mcarlo()

Dim i

Dim seed As LongLong

Dim x As LongLong
Dim period As Variant
Dim myarray() As LongLong
Dim counter As LongLong
Dim coord1 As LongLong
Dim coord2 As LongLong
Dim distance As LongLong
Dim pi As LongLong
Dim j As Variant
Dim k As Variant



period = 100
seed = 6


ReDim myarray(2 * period)

x = seed
For i = 1 To 2 * period
x = (4569872 * x + 1) Mod 45239875123

myarray(i) = x / 100000000

Next i

counter = 0
For j = 1 To period
coord1 = myarray(j)
For k = period + 1 To 2 * period
coord2 = myarray(k)
distance = Sqr(coord1 * coord1 + coord2 * coord2)
If distance <= 452.398 Then
counter = counter + 1
End If
Next k
Next j

pi = 4 * (counter / period)
MsgBox pi

End Sub

 

[StephenCrump]

Two observations:

I am not sure why you have a nested loop - for each random coord1, you're looping through the same 100 random coord2 values.
But if you're happy with the randomness of this approach, Pi will be 4 * counter / period ^ 2, rather than 4 * counter/period.

You also need to declare Pi as Double.

 

[KRice]

I agree with Stephen's observations. Several of the declarations need attention, as pi was being reported back as an integer. The nested loop is an issue, and I'm not sure about how the array values were being normalized...they should be divided by "m"...whichever value you choose in the mod function. Here's an edited version that gives pi=3.144 with different LCG parameters:

Option Explicit

Sub mcarlo()

Dim i
Dim seed As LongLong
Dim x As LongLong
Dim period As Variant
Dim myarray() As LongLong
Dim counter As LongLong
Dim coord1 As LongLong
Dim coord2 As LongLong
Dim distance As Variant
Dim pi As Double
Dim j As Long
Dim a As LongLong, c As LongLong, m As LongLong

period = 1000
'seed = 6

ReDim myarray(2 * period)
' LCG based on X_i = (a * X_(i-1) + c) mod m
a = 2175143
c = 10653
seed = 3553
m = 1000000

'Generate 2*period pseudo-random numbers
x = seed
For i = 1 To 2 * period
'x = (4569872 * x + 1) Mod 45239875123#
x = (a * x + c) Mod m

'myarray(i) = x / 100000000
myarray(i) = x

Next i


counter = 0
For j = 1 To period
coord1 = myarray(j)
coord2 = myarray(j + period)
'distance = Sqr(coord1 * coord1 + coord2 * coord2)
'If distance <= 452.398 Then
distance = Sqr((coord1 / m) * (coord1 / m) + (coord2 / m) * (coord2 / m))
If distance <= 1 Then
counter = counter + 1
End If
Next j

pi = 4 * (counter / period)
MsgBox pi

End Sub

 

[KRice]

Substituting back in your a,c,m parameters for the LCG, this gives pi~3.141756 with period=6,000,000 and seed=6. By normalizing the coord1 and coord2 with "m", then a distance <=1 lies inside the unit circle and would be picked up by the counter. Distances >1 lie outside the circle.

Option Explicit

Sub mcarlo()

Dim i
Dim seed As LongLong
Dim x As LongLong
Dim period As Variant
Dim myarray() As LongLong
Dim counter As LongLong
Dim coord1 As LongLong
Dim coord2 As LongLong
Dim distance As Double
Dim pi As Double
Dim j As Long
Dim a As LongLong, c As LongLong, m As LongLong

period = 6000000

ReDim myarray(2 * period)
' LCG based on X_i = (a * X_(i-1) + c) mod m
'a = 2175143
'c = 10653
'seed = 3553
'm = 1000000
a = 4569872
c = 1
m = 45239875123#
seed = 6

'Generate 2*period pseudo-random numbers
x = seed
For i = 1 To 2 * period
'x = (4569872 * x + 1) Mod 45239875123#
x = (a * x + c) Mod m

'myarray(i) = x / 100000000
myarray(i) = x
Next i

counter = 0
For j = 1 To period
coord1 = myarray(j)
coord2 = myarray(j + period)
'distance = Sqr(coord1 * coord1 + coord2 * coord2)
'If distance <= 452.398 Then
distance = Sqr((coord1 / m) * (coord1 / m) + (coord2 / m) * (coord2 / m))
If distance <= 1 Then
counter = counter + 1
End If
Next j

pi = 4 * (counter / period)
MsgBox pi

End Sub

 

[kisslaszlo123]

Thank you for your help Stephen and KRice. Now looking at the edited code, I understand what was the mistake. Thanks for the explanation too.