import math import numpy as np import random def w(x, mev, dis): return 1. / math.sqrt(2 * math.pi * dis) * math.exp(-(x - mev)**2 / (2 * dis**2)) def m(x): s = 0. for i in x: s += i return 1. / 20. * (s) def d(x, mev): s = 0. for i in x: s += (mev - i) ** 2 return math.sqrt(0.2 * s) x = [] for i in range(20): 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): x1 = [] maxW = w(mev, mev, dis) while 1: x = a + (b - a) * random.random() y = maxW * random.random() if y <= w(x, mev, dis): x1.append(x) if len(x1) == 20: return x1 x1 = Neumann(0, 6, mev, dis) mev1 = m(x1) dis1 = d(x1, mev1) print 'M-Neumann =', mev1 print 'D-Neumann =', dis1 def Muller(mev, dis): 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) == 20: return x2 print Muller(mev, dis)