pylab inline import math import numpy as np import matplotlib mlab as

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
%pylab inline
import math
import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import random
nsize = 20
def w(x, mev, dis):
return 1. / ( dis * math.sqrt(2 * math.pi)) * math.exp(-(x - mev)**2 / (2 * dis**2))
def m(x):
s = 0.
for i in x:
s += i
return 1. / nsize * (s)
def d(x, mev):
s = 0.
for i in x:
s += (mev - i) ** 2
return math.sqrt(1. / nsize * s)
x = []
for i in range(nsize):
x.append(np.random.random(6).sum())
mev = m(x)
dis = d(x, mev)
print 'M =', mev
print 'D =', dis
def Neumann(a, b, mev, dis, nsize):
x1 = []
maxW = w(mev, mev, dis)
while 1:
x = 6 * random.random()
y = maxW * random.random()
if y <= w(x, mev, dis):
x1.append(x)
if len(x1) == nsize:
return x1
x1 = Neumann(0, 6, mev, dis, nsize)
mev1 = m(x1)
dis1 = d(x1, mev1)
print 'M-Neumann =', mev1
print 'D-Neumann =', dis1
plt.hist(x1, 20)
plt.show()
m1 = []
d1 = []
for i in range(20):
x1 = Neumann(0, 6, mev, dis, nsize)
mmm = m(x1)
m1.append(mmm)
d1.append(d(x1, mmm))
mm1 = m(m1)
md1 = m(d1)
print mm1
print md1
print d(m1, mm1)
print d(d1, md1)
def Muller(mev, dis, nsize):
x2 = []
while 1:
x = random.uniform(-1, 1)
y = random.uniform(-1, 1)
s = x**2 + y**2
if 0 < s and s <= 1:
z0 = x * math.sqrt(-2 * math.log(s) / s)
z1 = y * math.sqrt(-2 * math.log(s) / s)
x2.append(mev + dis * z0)
x2.append(mev + dis * z1)
if len(x2) == nsize:
return x2
x2 = Muller(mev, dis, nsize)
mev2 = m(x2)
dis2 = d(x2, mev2)
print 'M-Muller =', mev2
print 'D-Muller =', dis2
plt.hist(x2, 20)
plt.show()
m2 = []
d2 = []
for i in range(20):
x2 = Muller(mev, dis, nsize)
mmm = m(x2)
m2.append(mmm)
d2.append(d(x2, mmm))
mm2 = m(m2)
md2 = m(d2)
print mm2
print md2
print d(m2, mm2)
print d(d2, md2)