B832@仲英楼 or 3133@泓理楼 // School of Physcis @ Xian Jiaotong University
一般来说，涉及到随机性的模拟方法都可称之为蒙特卡罗方法。通常也称为统计模拟方法。
Monte Carlo方法的名字来源于摩纳哥的一个城市蒙特卡罗，该城市以赌博业闻名。
MC 方法的提出可追溯到19世纪末期。20世纪40年代以后，随着电子计算机的发展，该方法得到迅速的发展和应用。
特别是在原子弹的研制过程中，Fermi, von Neumann, Ulam, Metropolis 等人为了研究核物理中的问题，如中子扩散等，系统发展了蒙特卡罗方法。
四位奠基人
除了在物理学及相关自然科学中的应用外，蒙特卡罗方法广泛的应用于社会学、金融等诸多领域。一个典型的例子是数量金融学。
随机性的例子 — 交通流模型
构造一个将积分区域包围起来的矩形 \([a, b] \times [0, y_0]\)。然后 随机的在矩形内撒 \(N\) 个点，假设其中有 \(M\) 个点落在积分线的下方，则积分为 \[ \int_a^b f(x)dx = \frac{M}{N}[y_0 (b - a)] \]
fx = lambda x: np.sin(2.0 * x)
x = np.linspace(0, 1, 100)
y = fx(x)
fig = plt.figure(figsize = (6, 3))
plt.plot(x, y)
N = 500000
xx = np.random.random(N)
yy = np.random.random(N)
yy_f = fx(xx)
indexes = np.where(yy > yy_f)
indexes2 = np.where(yy <= yy_f)
plt.plot(xx[indexes], yy[indexes], 'ro', markersize = 1)
plt.plot(xx[indexes2], yy[indexes2], 'bo', markersize = 1)
plt.show()
在区间 \([a, b]\) 内随机的撒 \(N\) 个点, 则 \(f(x)\) 在该区间内的积分可写为 \[ \int_a^b f(x)dx = \frac{(b - a)}{N}\sum_{i = 1}^N f(x_i) \qquad(1)\]
nn = 10000
M = 2**8
a = 57
c = 1
r1 = 10
def checkrand(r1, a, c, M, nn):
result1 = linmod(r1, a, c, M, nn)
result1 = result1 / M
plt.figure(figsize=(12, 3))
plt.subplot(131)
plt.plot(result1, '*', markersize = 1)
plt.xlabel(r'$i$')
plt.ylabel(r'$r_i$')
plt.subplot(132)
plt.plot(result1[0::2], result1[1::2], '*', markersize = 1)
plt.xlabel(r'$r_i$')
plt.ylabel(r'$r_{i+1}$')
plt.subplot(133)
c, xx = np.histogram(result1, np.linspace(0, 1, 50))
plt.plot(xx[1:] - 0.5 * xx[0], c, '-*')
plt.xlabel(r'$r_i$')
plt.ylabel(r'$N(r_i)$')
plt.suptitle('M = '+ str(M))
plt.show()
checkrand(r1, a, c, M, nn)
def checkrandpy():
result1 = np.random.random(10000)
plt.figure(figsize=(12, 3))
plt.subplot(131)
plt.plot(result1, '*', markersize = 1)
plt.xlabel(r'$i$')
plt.ylabel(r'$r_i$')
plt.subplot(132)
plt.plot(result1[0::2], result1[1::2], '*', markersize = 1)
plt.xlabel(r'$r_i$')
plt.ylabel(r'$r_{i+1}$')
plt.subplot(133)
c, xx = np.histogram(result1, np.linspace(0, 1, 50))
plt.plot(xx[1:] - 0.5 * xx[0], c, '-*')
plt.xlabel(r'$r_i$')
plt.ylabel(r'$N(r_i)$')
plt.suptitle('random from python')
plt.show()
checkrandpy()
N=np.array([10,20,50,100,200,500,1000,2000,5000,10000,20000,50000])
y = np.zeros(len(N))
err = y.copy()
sigma = y.copy()
def f(x):
return 1.0 / (1.0 + x**2.0)
for i in range(len(N)):
x = np.random.random(N[i])
fi = f(x)
favg = np.mean(fi)
y[i] = favg
err[i]=np.abs(y[i] - 0.25 * np.pi)
sigma[i] = np.sqrt(np.mean(fi**2.0) - favg**2.0) / np.sqrt(N[i])
plt.plot(N, err, 'r-x', label = 'error')
plt.plot(N, sigma, 'b-v', label = r'$sigma$')
plt.xscale('log')
plt.xlabel(r'$N_i$')
plt.title(r'error $\propto \sigma$ ')
plt.legend()
plt.show()
对前例 \(\int_0^1 \frac{dx}{1+x^2} = \frac{\pi}{4} \approx 0.78540\)
选择归一化的权函数 \(w(x) = \frac{1}{3}(4-2x)\)
其在积分区间上是单调下降的( \(f\) 也是)。由于在 \(x=0\) 和 \(x=1\) 两个端点上都有 \(f / w = 3 / 4\), \(w\) 对 \(f\) 的行为近似得很好。这样，新的积分变量为
\[ \begin{aligned} &y(x) = \int_0^x w(x')dx' = \frac{1}{3}x(4-x) \\ &\xrightarrow{\text{inverse function}} x = 2-(4-3y)^{1/2} \end{aligned} \]
\[ w(x) = \frac{6}{5}\left( 1 - \frac{1}{2}x^2\right) \rightarrow ~~~ y\sim x^3 \]
Ns = np.array([10,20,50,100,200,500,1000,2000,5000,10000,20000,50000])
res = np.zeros(Ns.shape)
for i in range(len(Ns)):
[x, y, z] = np.random.random((3, Ns[i]))
ff = np.exp(-0.5 * (x**2.0 + y**2.0 + z**2.0))
res[i] = (2.0 * np.pi)**(-1.5) * np.mean(ff)
plt.plot(Ns, res - 0.039772182)
plt.xscale('log')
plt.show()
b = 1.0
a = 0.0
w = lambda x: 1.2 * (1.0 - 0.5 * x * x)
ymax = np.max(w(np.linspace(a, b, 1000)))
nn = 1000000
i = 0
xres = np.zeros(nn)
while i < nn:
xx = a + (b - a) * np.random.random()
yy = ymax * np.random.random()
if (yy < w(xx)):
xres[i] = xx
i += 1
cc, cx = np.histogram(xres, bins = 100)
dx = cx[1] - cx[0]
plt.plot(cx[1:] - 0.5 * (cx[1] - cx[0]), cc / nn / dx)
plt.plot(np.linspace(a, b, 100), w(np.linspace(a, b, 100)))
plt.show()
def metroplis1d(x0, delta, w):
xt = x0 + delta * (2.0 * np.random.random() - 1.0)
wx0 = w(x0)
wxt = w(xt)
if (wxt > wx0):
return xt
else:
prob = np.random.random()
if (wxt / wx0 > prob):
return xt
else:
return x0
x0 = 0
nn = 10000
xres = np.zeros(nn)
i = 0
xres[0] = 0
wfunc = lambda x: 1.0 / np.sqrt(2.0) * np.exp(-0.5 * x * x)
while i < nn-1:
xres[i+1] = metroplis1d(xres[i], 2.0, wfunc)
i += 1
plt.hist(xres, bins = 100)
plt.show()
def metroplis2d(x0, delta, w):
xt = x0 + delta * (2.0 * np.random.random(2) - 1.0)
wx0 = w(x0)
wxt = w(xt)
if (wxt > wx0):
return xt
else:
prob = np.random.random()
if (wxt / wx0 > prob):
return xt
else:
return x0
x0 = 0
nn = 100000
xres = np.zeros([nn, 2])
i = 0
xres[0] = 0
wfunc = lambda x: 1.0 / np.sqrt(2.0) * np.exp(-0.5 * (x[0] * x[0] + x[1] * x[1]))
while i < nn-1:
xres[i+1, :] = metroplis2d(xres[i, :], 2.0, wfunc)
i += 1
plt.hist2d(xres[:, 0], xres[:, 1], bins = 100);
plt.show()
def metroplis1d(x0, delta, w):
xt = x0 + np.random.random() - 0.5
wx0 = w(x0)
wxt = w(xt)
if (wxt > wx0):
return xt
else:
prob = np.random.random()
if (wxt / wx0 > prob):
return xt
else:
return x0
x0 = 0.5
nn = 1000000
xres = np.zeros(nn)
i = 0
xres[0] = 0
wfunc = lambda x: 1.2 * (1 - 0.5 * x * x) if (x <= 1 and x >= 0) else 0.0
wfuncc = lambda x: 1.2 * (1 - 0.5 * x * x)
while i < nn-1:
xres[i+1] = metroplis1d(xres[i], 0.5, wfunc)
i += 1
cc, cx = np.histogram(xres, bins = 100)
plt.plot(cx[1:] - 0.5 * (cx[1] - cx[0]), cc / (cx[1] - cx[0]) / nn)
plt.plot(np.linspace(0, 1, 100), wfuncc(np.linspace(0, 1, 100)))
plt.show()
