Sampling Functions for PI

'''
Functions From Statistical Mechanics 
by Werner Krauth

'''

import random
import math

def directsampling(n,r):
    data=[]
    print 'Direct Sampling Case'
    print 'Runs       Number Hits  estimate of pi'
    for jj in range(r):
        nhits = 0
        for ii in range(n):
            x=random.uniform(-1,1)
            y=random.uniform(-1,1)
            if (x**2+y**2) < 1:
                nhits= nhits + 1
        print str(jj)+'            '+str(nhits)+'               '+str(4.0*nhits/n)
        data.append([jj,nhits])
    return data

def markovPI(n,r,throwdist):
    data=[]
    print 'Markov Chain Case distance = ' +str(throwdist)
    print 'Runs       Number Hits  estimate of pi'
    for jj in range(r):
        nhits = 0
        x = 1
        y = 1
        for ii in range(n):
            dx = random.uniform(-throwdist,throwdist)
            dy = random.uniform(-throwdist,throwdist)
            if (abs(x+dx) < 1 and abs(y+dy)<1):
                x=x+dx
                y=y+dy
            if (x**2 + y**2) < 1:
                nhits= nhits + 1
        print str(jj)+'            '+str(nhits)+'               '+str(4.0*nhits/n)
        data.append([jj,nhits])
    return data

def buffon(n,r,a,b):
    '''
	Buffon's Needle Experiment
	n = number of throws
	r = number of runs
	a = length of needle
	b = distance between cracks
	theta = angle needle makes to crack
	rcenter = center of needles on floor
	0  < theta < pi/2
	0 < xcenter < b/2

nhits <===  number of hits of needle centered at x, with orientation theta
	nhits = 1 if x < a/2 and abs(theta) < arcos(x/(a/2))
      = 0 otherwise
	'''
    data=[]
    print 'Buffon Needle Experiment (Google it) ' 
    print 'Runs       Number Hits  estimate of pi'
    for jj in range(r):
        nhits = 0
        for ii in range(n):
            xcent = random.uniform(0,b/2.0)
            theta = random.uniform(0,math.pi/2)
            xtip  = xcent - (a/2.0)*math.cos(theta)  #use of cosine not historically accurate
            if xtip < 0 :
                nhits += 1
        #print str(jj)+'            '+str(nhits)+'               '+str((6.0/a*float(b))*nhits/n)
        c = 2.0*a*n
        d = b*nhits
        print str(jj)+'            '+str(nhits)+'               '+str(c/d)
        data.append([jj,nhits])
    return data
    
def buffon1(n,r,a,b):
    '''
    Buffon's Needle Experiment
	n = number of throws
	r = number of runs
	a = length of needle
	b = distance between cracks

nhits <===  number of hits of needle centered at x, with orientation theta
    This buffon's need is more historically accurate.  It does not use pi nor does 
    it use the cosine function
	'''
    data=[]
    print 'Buffon Needle Experiment (Historical) ' 
    print 'needle = '+str(a) + '  crack = '+str(b)
    print 'Runs       Number Hits   est. of PI    '
    for jj in range(r):
        nhits = 0
        for ii in range(n):
            xcent = random.uniform(0,b/2.0)
            tau=2
            while tau > 1 :
            	dx = random.uniform(0,1)
            	dy = random.uniform(0,1)
            	tau = math.sqrt(dx**2 + dy**2)
            	if tau <= 1:
                    xtip  = xcent - (a/2.0)*dx/tau
                    if xtip < 0 :
                        nhits += 1
       
        c = 2.0*a*n
        d = b*nhits
        print str(jj)+'             '+str(nhits)+'                '+str(c/d)
        data.append([jj,nhits])
    return data

r=5
n=1000
hits=directsampling(n,r)
dist = 2  #adjust this and observe variations
hits=markovPI(n,r,dist)
a = 2  #needle 2 inches
b = 3  #cracks 1 inch spacing
hits= buffon(n,r,a,b)
hits= buffon1(n,r,a,b)

Python Fiddle

Python Cloud IDE