Algorithm 1.1 Buffon

throws = 10000
x = np.random.rand(throws) # 10000 vector of samples from U(0,1)
theta = 0.5*np.pi*np.random.rand(throws) # 0 < theta's < pi/2
hits = (x<0.5*np.sin(theta)) # vector of true/false
print 'pi estimate:',throws/float(sum(hits))

Algorithm 1.2 Buffon Ab Initio

# Ab initio buffon needle program to estimate pi by Buffon simulation
# without knowing pi to start with.
# To convert to the simulation with pi known, use PI instead of pi.

import numpy as np
trials = 10000
segments = 20 #avoid division by 0

PI = np.pi
d = 2; L = 1; #for reference
pi = 1.0 # starting value for iteration
throws = 0; hits = 0

for i in range(trials):
  x = np.random.rand(segments) # 0 < x < d/2
  theta = 0.5*pi*np.random.rand(segments) # 0 < theta < pi/2
  hits += (x<0.5*np.sin(theta)) # unlikely to be 0 after 20 throws
  throws += segments
  pi = throws/float(sum(hits))
print('estimate of pi after', throws, 'throws:', pi)

Algorithm 1.3 Buffon's Histogramming Algorithm

import numpy as np
import matplotlib.pyplot as plt
throws = 100 # 100 trials
replications=4000
y= [0.]*replications
for i in range(replications):
 x=np.random.rand(throws)
 theta=0.5*np.pi*np.random.rand(throws)
 hits= (x<=0.5*np.sin(theta))
 y[i] = (sum(hits)/float(throws)) # y= result vector
recipEst=sum(y)/float(replications)
vVec=(y-recipEst)*(y-recipEst); # vector of sqr deviations
print vVec
v=sum(vVec)/(replications-1); # sample variance
stddev=np.sqrt(v) # sample standard deviation
piEst=1/recipEst
print piEst
#fig = plt.figure(figsize=(6,6))
fig = plt.figure()  # defaults to 640x480, 100 px/in
plt.hist(y,bins=12)
#fig.savefig('histoPlot.png', bbox_inches='tight', dpi=150)
fig.savefig('recipHist.png')
plt.show()

Algorithm 1.4 Gambler's Fortune

import numpy as np
import random
import matplotlib.pyplot as plt

R = []
i=0; R.append(100); # gambler's initial fortune}
while( R[i] > 0 and R[i]<2100 ):
  i=i+1
  W = (random.random() < 0.5); # random value of 0 or 1}
  W = 2*W - 1; # random value of +/-1}
  R.append(R[i-1]+W) # gamblers new fortune}

x = np.linspace(1,len(R),len(R))
plt.plot(x,R) # plot R against its index
plt.show()
# run for investigating running time
T = []
for j in range(20):
  R = []
  i=0; R.append(100);
  while( R[i] > 0 and R[i] < 2100 ):
    i=i+1
    W = 2*(random.random() < 0.5)-1
    R.append(R[i-1]+W)
  T.append(i)
fig = plt.figure(figsize=(6,6))
#fig = plt.figure()  # defaults to 640x480, 100 px/in
plt.hist(T,bins=12)
#fig.savefig('histoPlot.png', bbox_inches='tight', dpi=150)
#fig.savefig('histoPlot.png')
plt.show()

Algorithm 1.5 Dice Rolling

import numpy as np
import random
import matplotlib.pyplot as plt

nRolls = 60
red = np.array([random.randint(1,6) for i in range(nRolls)])
green = np.array([random.randint(1,6) for i in range(nRolls)])
dice = red+green
b=round(1+13/0.1)
x = np.linspace(0,13,b)
cdf = []
for i in range(len(x)):
  cdf.append(sum(dice <= x[i]))
fig = plt.figure()
plt.title("Dice pair cdf for "+str(nRolls)+" rolls")
plt.plot(x,cdf)
fig.savefig('dicerollCDF.png')
plt.show()

Algorithm 1.6 Joint Probabity Histogram

import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

mat = np.random.randint(low = 10, high = 30, size = (60,2))
#print(mat)
df4 = pd.DataFrame(mat,columns=["x","y"])

sns.set(font_scale=1.44) # increase from standard 1
sns.jointplot(data=df4,x="x",y="y",kind='hist',bins=20,
  marginal_kws=dict(bins=20))
## each blue square = 1/60, black = 2/60
#sns.jointplot(data=df4,x="x",y="y",kind='hist',marginal_kws=dict(fill=False))
#  marginal_kws=dict(bins=25, fill=False))

plt.suptitle("Joint Histogram and Marginals", y=1.01) # y raises the title
#plt.savefig("jointHistogram.png")
plt.show()

Algorithm 1.7 Middle Sine

import numpy as np
import math
import time
import random

tenpow5 = 100000
nIts = 1000000
x= 0.5
for i in range(10):
  x = tenpow5*np.sin(x)
  w = math.floor(x)
  x = x-w
  print(x)
# timing
start = time.time()
for i in range(nIts):
  x = tenpow5*np.sin(x)
  w = math.floor(x)
  x = x-w
end = time.time()
print('time for '+str(nIts)+' middle digits',end-start) # 2.03615999222

start = time.time()
for i in range(nIts):
  #x = np.random.rand(1); # 0.814379930496
  x = np.random.random_sample(); # 0.28878903389
end = time.time()
print('time for '+str(nIts)+' built-in in a loop ',end-start)

start = time.time()
y = np.random.rand(nIts) # 0.0099949836731
end = time.time()
print('time for '+str(nIts)+' vector generation',end-start)

Algorithm 1.8 Middle Square

import math
tenpow7 = 10000000
x = 0.1234567; # initial seed
for i in range(2):
   x = x*tenpow7; # x now an integer
   w = x*x;  # square
   w = math.floor(w/1000); # drop last three digits
   w = w/tenpow7;  # 7 decimal digits
   x = w - math.floor(w) # drop the integer digits
   print('next seed ',str(x))

Algorithm 1.9 Hit-or-Miss for pi

import numpy as np
nTrials = 10000
x = np.random.rand(nTrials)
y = np.random.rand(nTrials)
hits=sum(x**2+y**2<1) # x squared plus y squared
piEst = 4*hits/float(nTrials)
print piEst

Algorithm 1.10 Coupon Collecting

import numpy as np
import random
import matplotlib.pyplot as plt

nLetters = 9  # ANGELFISH}
nTrials = 10000
nTries = [0]*nTrials
for i in range(nTrials):
  success = 0
  nTries[i] = 0  # redundant but reassuring
  angelfish = [0]*nLetters # reset no letters achieved
  while success == 0:
    nTries[i] += 1  #inc. count
    buy = random.randint(0,nLetters-1) #letter obtained
    angelfish[buy] = 1
    if sum(angelfish)==nLetters:
      success = 1
#print nTries
fig = plt.figure()
plt.hist(nTries)
plt.show()

Algorithm 1.11 Permutations

import random
n = 9 # a permutation of 9 objects
m = 5 # do 5 of them
perm = [i for i in range(n)]  #identity permutation
for i in range(m):
  j = random.randint(i,n-1) # from the remaining
  swap = perm[i]
  perm[i] = perm[j]
  perm[j] = swap
  result = [perm[i] for i in range(m)]
print result