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) ```