CODEWITHSUNDEEP
pythonmathpi
nuwanda123
nuwanda123

Reputation: 17

Estimating pi with a Monte Carlo method results in a larger value than expected

I'm trying to estimate pi by dividing the areas of a square and its embedded circle, but I get ~3.66.

Does anyone see what I'm doing wrong?

inCount=0
outCount=0
it=1000000
L=100
for i in range(it):
    xran=rnd.random()*L
    yran=rnd.random()*L
    xc=abs(0.5*L-xran)
    yc=abs(0.5*L-yran)
    r=np.sqrt((xc**2)+(yc**2))
    if r<0.5*L:
        inCount=inCount+1
    if r>0.5*L:
        outCount=outCount+1
    if r==0.5*L:
        inCount=inCount+1
        outCount=outCount+1
pigen=inCount/outCount
print('pi generated: '+np.str(pigen))

Upvotes: 1

Views: 159

Answers (4)

Matthias Fripp
Matthias Fripp

Reputation: 18625

inCount shows the number of points within a circle of radius r = L/2 and outCount shows the number of points in a square exactly containing that circle, but not in the circle itself.

inCount is proportional to pi * r**2 and outCount is proportional to L**2 - pi * r**2 = (4 - pi) * r**2. When you take the ratio you get pi / (4 - pi) = 3.66.

As Wai Ha Lee pointed out, you need to calculate 4 * inCount / (inCount+outCount) = pi.

Upvotes: 0

Wai Ha Lee
Wai Ha Lee

Reputation: 8805

You have

pigen=inCount/outCount

which gives the proportion of hits inside to outside the radius.

Note that pi/(4-pi) = 3.659792... which is what your code is currently estimating.

You need

pigen=4*inCount/(inCount+outCount)

which will give you four times the proportion of hits inside compared to the total, i.e. pi.

Also note that your code is currently

if r<0.5*L:
    inCount=inCount+1
if r>0.5*L:
    outCount=outCount+1
if r==0.5*L:
    inCount=inCount+1
    outCount=outCount+1

which can be simplified with elif/else. Since r cannot be both greater than, and less than, L, the second if can become an elif. Likewise, if r is neither less than, or greater than, L then it must be equal so the third if can simply become an else.

if r<0.5*L:
    inCount=inCount+1
elif r>0.5*L:
    outCount=outCount+1
else:
    inCount=inCount+1
    outCount=outCount+1

This will prevent unnecesary comparisons of r and L in your code.

Your final code would then be

inCount=0
outCount=0
it=1000000
L=100
for i in range(it):
    xran=rnd.random()*L
    yran=rnd.random()*L
    xc=abs(0.5*L-xran)
    yc=abs(0.5*L-yran)
    r=np.sqrt((xc**2)+(yc**2))
    if r<0.5*L:
        inCount=inCount+1
    elif r>0.5*L:
        outCount=outCount+1
    else:
        inCount=inCount+1
        outCount=outCount+1
pigen=pigen=4*inCount/(inCount+outCount)
print('pi generated: '+np.str(pigen))

Upvotes: 7

nuwanda123
nuwanda123

Reputation: 17

Wai Ha Lee was almost right! I also had forgotten to add a 4.

If anyone wants to know how it goes, it's like that (and yes, L can be whatever you want):

import numpy as np
import random as rnd
inCount=0
outCount=0
it=1000000
L=100
for i in range(it):
    xran=rnd.random()*L
    yran=rnd.random()*L
    xc=abs(0.5*L-xran)
    yc=abs(0.5*L-yran)
    r=np.sqrt((xc**2)+(yc**2))
    if r<0.5*L:
        inCount=inCount+1
    if r>0.5*L:
        outCount=outCount+1
    if r==0.5*L:
        inCount=inCount+1
        outCount=outCount+1
pigen=4*inCount/(inCount+outCount)
print('pi generat: '+np.str(pigen))

Upvotes: 0

Jiaju Huang
Jiaju Huang

Reputation: 79

inCount+outCount = 4*r^2
inCount = pi*r^2

So if you need to get pi

pigen=inCount/(outCount+inCount)*4

Upvotes: 1

Related Questions