Monte Carlo Integration
-
Approximating the area under a standard normal curve.
import scipy from scipy.stats import distributions as ssd def test3_1(): m = 10000 a,b = 0., 1. w = (b-a)/m x = np.linspace(a+w, b-w, m) y = np.exp(-x*x/2)/np.sqrt(2*np.pi) print (y*w).sum() def test3_2(): m = 10000 a,b = 0., 1. w = (b-a)/m x = (b-a)*np.random.random_sample(size=m)+a y = ssd.norm.pdf(x) print (y*w).sum() -
Acceptance-Rejection Method
def test3_3(): sz = 500000 u = np.random.random_sample(size=sz) h = np.random.random_sample(size=sz)*0.4 # point (u,h) num = (h<ssd.norm.pdf(u)).mean() print num*0.4 -
Random Sampling Method
def test3_4(): m = 500000 a, b = 0, 1 y = ssd.norm.rvs(size=m) print ((a<=y)&(y<=b) ).mean()
Box-Muller method
This method and variations of it are widely used in statistical software to simulate random observations from NORM(0; 1)
def test2_6():
m,bins = 10000, 50
u1 = np.random.random_sample(size=m)
u2 = np.random.random_sample(size=m)
a = np.sqrt(-2*np.log(u1))
b = 2*np.pi*u2
z1 = a*np.cos(b)
z2 = a*np.sin(b)
pl.subplot(131)
pl.hist(z1, bins=bins)
pl.subplot(132)
pl.hist(z2, bins=bins)
pl.subplot(133)
pl.ylim(0,1)
hy,hx = pl.histogram(z1,bins=bins)
pl.bar(hx[:-1],hy*1.0/max(hy)/np.sqrt(2*np.pi), width=1.0/(len(hx)-1),color='green' )
exponential distribution
simulate an observation X that has an exponential distribution.
def test2_5():
m,bins = 10000, 50
u = np.random.random_sample(size=m)
x = -np.log(1-u)
pl.subplot(121)
pl.hist(x, bins=bins)
pl.subplot(122)
pl.plot(x,np.exp(-x), '.', color='red' )
Square of a uniform random variable
def test2_4():
m,bins = 100000, 50
u = np.random.random_sample(size=m)
x = u**2
alpha = 1.0*m/bins
pl.subplot(121)
pl.ylim(0,alpha+100)
pl.hist(u, bins=bins)
pl.subplot(122)
pl.ylim(0,8)
hy,hx = pl.histogram(x,bins=bins)
pl.bar(hx[:-1],hy/alpha, width=1.0/(len(hx)-1),color='green' )
pl.plot(x,0.5*x**(-0.5), '.', color='red' )
.png)
UNIF(0,1) generator
def test2_3():
a,b,d = 1093,18257, 86436
s, m = 7, 1000
rs = np.zeros(m, int)
rs[0] = s
for i in xrange(m-1):
rs[i+1] = (rs[i]*a+b)%d
u = (rs+0.5)/d
pl.subplot(121)
pl.hist(u)
pl.subplot(122)
pl.plot(u[:m-1],u[1:m], '.')
Birthday Matches Combinatorial Approach
Suppose there are n =25 people in a room.What is the probability that two or more of them have the same birthday?
#coding=utf-8
import numpy as np
import pylab as pl
import time
np.random.seed(int(time.time()))
def birthday():
m = 60
x = xrange(1,m)
y = list(1-np.prod([1-i/365. for i in range(e)]) for e in x)
pl.plot(x,y,'-')
if __name__=='__main__':
birthday()
pl.show()
Sampling Computer Chips
Suppose there are 100 memory chips in a box, of which 90 are good and 10 are bad. We withdraw five of the 100 chips at random to up grade a computer.What is the probability that all five chips are good?
#coding=utf-8
import numpy as np
import time
np.random.seed(int(time.time()))
def Sampling():
num, pkg, sz, iGood = 10000,100,5,90
good= np.zeros(num,int)
x= np.array(xrange(pkg))
for i in xrange(num):
np.random.shuffle(x)
e = x[:sz]
count= (e<iGood).sum()
good[i] = count
print (good==sz).mean()
answer: 
