Algorithm 2.1 Invert Multi-outcome Discrete Density

N = 4 # number of outcomes
oc = [2,3,4,5] # list the outcomes, 2 through 5 eggs
op = [0.15,0.2,0.6,0.05] # list of outcome probabilities
bp = [0]*N # set up the breakpoints
bp[0]=op[0]
for i in range(1,N):
   bp[i] = bp[i-1] + op[i]
# end of set up, go here to get samples
U=random.random(); k=0;
while( U >= bp[k]):
   k=k+1;
selected = oc[k]
print(selected)

Algorithm 2.2 Roulette Wheel Selection

import random
N = 4
tabLen = 100
oc = [2,3,4,5]
op = [0.15,0.2,0.6,0.05] # event probabilities
op =[int(op[i]*tabLen) for i in range(len(op))] # 15,20,60,5
 # set up: assign fractions of 100 according to the probabilities
rv = []  # for return values
for i in range(len(op)):
    for j in range(op[i]):
      rv.append(oc[i])
 # end of setup, go here to get samples
k = random.randint(0,99) # rv[k] is selected
selected = rv[k]

Algorithm 2.3 Vector Algorithm for Discrete Selection

import random
op = [0.15,0.2,0.6,0.05]}
bp = []
bp.append(op[0])
for i in range(1,len(op)):
  bp.append(bp[i-1]+op[i])
N = 10000
u = np.random.rand(N) # N vector of samples from U(0,1)
w2 = (u < bp[0])*2  # 2's where true
a = np.array(u>=bp[0]) # working arrays needed
b = np.array(u< bp[1])  # for the and operation
w3 = (a & b)*3  # 3's where true
a = np.array(u>=bp[1])
b = np.array(u< bp[2])
w4 = (a & b)*4  # 4's where true
w5 = (u>=bp[2])*5 # 5's where true
w = w2+w3+w4+w5 # 10000 vector with the right frequencies

Algorithm 2.4 Walker Construction Algorithm

import numpy as np
import matplotlib.pyplot as plt

def firstUnder(arr): #but greater than 0
  k=0
  while( True ):
    if( arr[k] == 0 or arr[k] > S ):
      k += 1
    else:
      break;
  return k

def firstOver(arr):
  k = 0
  while( True ):
    if( arr[k] < S):
      k += 1
    else:
      break;
  return k

def reset(idx,ind):
  Q[idx] = tab[idx]
  tab[idx] = 0
  resid = S - Q[idx]
  A[idx] = ind
  tab[ind] -= resid

### mmm ###
#rw = [.3,.4,.1,.2]
rw = [.02,.04,0.08,.15,.2,.21,.16,.09,.05];
M = len(rw)
tab = [0]*M
for i in range(M):
  tab[i] = rw[i]

Q = [0]*M
A = [0]*M
S = 1.0/M

for i in range(M-1):
  t = firstUnder(tab) #take from this sector
  p = firstOver(tab)  #put the residual here
  reset(t,p)

#last
t = 0; max = 0;
for i in range(M):
  if( tab[i] > max):
    max = tab[i];
    t = i
Q[t] = tab[t]


#print("Q-pre",Q)
for i in range(len(Q)):
  Q[i] *= len(Q)
print("Q: ",end="")
for i in range(len(Q)):
  print(str(round(Q[i],3))+" ",end="")
print("")
print("A: ",A)

Algorithm 2.5 Poisson Sampling

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

def poissonSample(targ):
  u = 1; x = 0
  while( u > targ):
    last = x
    u *= random.random()
    #print("u now "+str(u)+", last now "+str(last))
    x += 1
  return last

lam = 0.1 # event rate per minute, lambda is a reserved word
t = 30  # number of minutes
targ = np.exp(-lam*t)
 # get a sample
X = poissonSample(targ)
print(X)

Algorithm 2.6 Rational Approximation to the Normal CDF

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

# Abramowitz and Stegun
def normalCDF(x):
  a1 = 0.31938153; a2 = -0.356563782;
  a3 = 1.781477937; a4 = -1.821255978;
  a5 = 1.330274429; b=0.2316410; c=1/np.sqrt(2*np.pi)
  L=abs(x); K=1/(1+b*L);
  w = 1-c*np.exp(-L*L/2)*K*(a1+K*(a2+K*(a3+K*(a4+a5*K))));
  msign = -(x<0); mfact = 2*msign+1;
  w = mfact*(w+msign)
  return w
x = 1
print(normalCDF(x))

Algorithm 2.7 Gilbrat's Approximation Example

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

a = 0.2; b=2.0; n=10
c = b-a
alfa = (n/c)*(b*(np.log(b)-1)-a*(np.log(a)-1))
beta2 = (n/c)*(b*(np.log(b)*np.log(b)-2*np.log(b)+2)-\
a*(np.log(a)*np.log(a)-2*np.log(a)+2))-alfa**2
beta = np.sqrt(beta2)

nReps = 2000000
X = np.zeros(nReps)
for i in range(nReps):
  R = 1;
  for j in range(n):
    r = a + c*random.random()
    R *= r
  X[i] = R

## generate the normal samples using alfa and beta
logXexact = beta*np.random.randn(nReps)+alfa

# do the lognormal
Xexact = np.exp(logXexact)
xmin = 0; xmax = 200

fig = plt.figure()
ax = plt.subplot(111)
ax.hist(Xexact, bins=200, density=True, log = True, \
range = [0,xmax],label="exact lognormal")
plt.hist(X, bins=200, density=True, log = True, \
range = [0,xmax],label="product random variable")
ax.legend()
x = np.linspace(xmin,xmax)
y = (1/(x*beta*np.sqrt(2*np.pi)))*\
np.exp(-(np.log(x)-alfa)**2/(2*beta*beta))
ax.plot(x, y, 'k', linewidth=1) # plot in black
plt.xlabel("X")
plt.ylabel("probability density")

title = "Lognormal Histogram Comparison, alpha= {:.2f}, beta= {:.2f}".\
format(alfa, beta)
plt.title(title)
plt.show()
fig.savefig('lognormal_comparisonDisp.png')

Algorithm 2.8 Joint Sampling Over a Triangle

mport numpy as np
import matplotlib.pyplot as plt

x = [0,0,1,0]
y = [0,1,0,0]
plt.plot(x,y)

for i in range(200):
  u1 = np.sqrt(np.random.uniform())
  u2 = np.random.uniform()
  x = 1-u1
  y = u2*u1
  plt.scatter(x,y,color='black',s=8)

plt.text(0.40,0.67, "y=1-x", fontsize=12)
plt.title("200 points uniformly over a triangle")
plt.hlines(0.685,0.322,0.39)
#plt.savefig("uniformovertrianglefig.png")
plt.show()