不会被那透明的风暴所迷惑，一定要找到你。——ユリ熊嵐
关于蒙特卡洛
这个学期上统计，学的时候提到蒙特卡洛方法。这个名字听起来很高级，其实就是上随机数乱搞(bushi)。对于一个不熟悉的分布，工程上可以通过生成符合该分布的随机数的方法研究其相关性质，在大数律的保证下，只要生成的随机数足够多，总是能达到需要的精度的。方法很朴实，效果却不错。至少对于像我这样的数学苦手而言，这个方法在很多时候可以避免对付复杂的积分。或者是在研究一些比较冷门的性质的时候，这个方法也可以用于直观地计算。
举个例子Problem #3658 - ECNU Online Judgeacm.ecnu.edu.cn
若干个月前的EOJ月赛题，大致上是求以下内容：
将
个点随机地投射到一个大小为
的矩形中，每个点的横纵坐标满足独立同分布
，求点两两之间最近距离的期望是多少。
这个题要是手算人就没了，但用蒙特卡洛就很好理解，生成均匀分布的点算最小距离即可。
话归正传，可以总结蒙特卡洛的基本方法：判断需要什么分布
生成符合分布的随机数
对随机数进行计算
除了计算这些随机数的特征以外，蒙特卡洛的另一个用途就是可以结合可视化，绘制一些分布的直观表示，比如分布的直方图等等。
听起来似乎很简单，具体实现起来也不会太复杂。下面以生成正态分布的动图为例，看一下怎么具体实现蒙特卡洛方法。
绘制动图
基于python实现很简单，因为我们有万能的plt。在这里高呼三声：
PLT NB！
PLT NB！
PLT NB！
然后开始绘制，直方图的绘制是很容易的，直接使用 plt.hist 就可以了。比较麻烦的是生成动图。一开始的想法是先生成足够多张图，然后再调库合并成gif。这么复杂的方法(还有需要中间输出)明显是绕了弯路，但或许是早上刚睡醒神志不清，我在看着程序跑的时候才感觉不对，打开搜索引擎，发现matploylib有提供合成动画的模块animation。对着代码实现
代码大概长这样子：
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
import numpy as np
fig = plt.figure()
data = np.array(0)
def update(frame):
plt.cla()
plt.ylim(0,100)
plt.xlim(-5,5)
global data
data = np.append(data , np.random.randn(20))
hist = plt.hist(data, 60,density = False,facecolor="blue")
return hist,
ani = FuncAnimation(fig, update, frames = 100,interval=200)
ani.save("./test-1.gif",writer='pillow')
plt.show()采用的是animation里面的FuncAnimation方法
这个动图是蒙特卡洛在可视化上的应用，如果想要利用蒙特卡洛计算统计量，当然也是可以的。为了同时体现蒙特卡洛的两种应用，在原来的代码上再加一点计算均值和方差的内容，同时表示在动图里。
可以看到，均值和方差与采用的分布(标准正态)基本吻合，可见生成的随机数确实可以用来计算统计量。虽然均值和方差是容易知道的，似乎没有使用蒙特卡洛的必要，但其他的统计量也是同样的道理。
具体的代码实现如下：
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
import numpy as np
fig = plt.figure()
data = np.array(0)
def update(frame):
# 清除画布，不然均值和标准差的线会糊成一团
plt.cla()
plt.ylim(0,100)
plt.xlim(-5,5)
# 不断追加生成的随机数
global data
data = np.append(data , np.random.randn(20))
# 计算均值
mean = np.mean(data)
plt.axvline(x=mean,ls="-",label="mean")
# 计算标准差，绘制3 sigma的范围
std = np.std(data)
plt.axvspan(xmin = mean - std,xmax = mean + std,facecolor = 'blue',alpha=0.2,label = "1 $\sigma$")
plt.axvspan(xmin = mean - 2*std,xmax = mean + 2*std,facecolor = 'blue',alpha=0.2,label = "2 $\sigma$")
plt.axvspan(xmin = mean - 3*std,xmax = mean + 3*std,facecolor = 'blue',alpha=0.1,label = "3 $\sigma$")
plt.legend()
# 绘制图像
hist = plt.hist(data, 60,density = False,facecolor="blue")
return hist,
ani = FuncAnimation(fig, update, frames = 100,interval=200)
ani.save("./test.gif",writer='pillow')
plt.show()生成各种分布的随机数
在上文中，对于怎么生成随机数一直是一笔代过。样例也是最简单的正态分布，任何讲python生成随机数的文章都必会涉及的内容，numpy也提供了简单的调库生成方法。那如果要实现比较复杂的分布，难道还能调库吗？
在大部分情况下，答案是能的，因为numpy提供了大量分布的随机数生成方法见numpy文档
如果还是不行，那也不是没有办法：如果已知密度函数，可以生成均匀分布的随机数，用变量变换搞一下；
如果不知道密度函数，但知道生成方法(e.g. 抽卡机制)，直接逐个模拟；
如果什么都不知道，那可以强行上渐进正态(并不能)，或者退出这篇文章，假装不存在这个方法。
总结
这里主要水了一下蒙特卡洛的实现，顺带介绍了用plt绘制动图的方法。
一句话概括蒙特卡洛，蒙特卡洛就是生成随机数乱搞。
一句话概括绘制动图，调用FuncAnimation，在update函数里写绘制代码。
