用蒙特卡洛求pi_python 实现蒙特卡洛方法

本篇来讲一个简单的蒙特卡洛方法及其应用

什么是蒙特卡洛方法

蒙特卡洛方法是一种在统计学十分常用的一种随机模拟方法, 简单来说即使用随机数来解决计算问题.

在了解蒙特卡洛方法之前, 我们先要回顾几个统计学概念:

1. 期望(平均数)
2. 方差, 标准差 (
   )
3. 正态分布

即一个随机变量x服从一个数学期望

、方差为

的概率密度函数可以呈现出图1中的钟形曲线. 期望

决定了曲线的位置(比如图中蓝色和红色曲线

); 标准差

决定了分布的幅度,

越大代表数据分布越广. 在这里我们还需要知道三西格玛定律, 即

也就是说,当标准差取

的时候,数据的概率密度已经可以覆盖近乎99%了.

图1, 引自wiki

回顾了这个概念后我们不难理解, 在我们日常处理的分析数据中, 我们利用n多个数据求得的平均数和标准差其实只是代表了这个数据的概率密度分布情况. 在这里我们举个简单的例子来应用蒙特卡洛方法. 我们用python来实现这个过程:
首先我们自定义一个变量x的平均值为10, 标准差为0.3

import numpy as np
import random
import matplotlib.pyplot as plt

meanX, stdX=10,0.3

现在我们打算根据正态分布的公式来随机产生1000次x可能出现的数值

n=1000 #这里的n就是蒙特卡洛模拟的随机数生成器

x=np.random.normal(meanX, stdX, n) #使用numpy内置的正态分布函数random.normal(),随机产生x1000次

我们再假设自变量x和因变量y有如下函数关系:

y=x*10+x**2+5
np.shape(y) #用numpy的shape函数可以查看变量的大小,结果显示y为一个1000行的数组
(1000,)

计算取得y的平均值和标准差

meanY=np.mean(y)
stdY=np.std(y)
meanY
205.26338508878402
stdY
9.301438243719513

由此我们得出了变量y在随机模拟1000次下可能产生的数值,平均值大约为204.9,而标准差已经扩大到9.2了.换句话说y的概率分布也更大了,这也意味着由y再参与的计算已经不像x那样具有更收敛的计算结果了.

# 查看一下x和y随机产生了那些数值
fig,axes=plt.subplots(1,2,figsize=(8,4),constrained_layout=True)
axes[0].bar(np.arange(n),x)
axes[1].bar(np.arange(n),y)
axes[0].set_xlabel('n')
axes[0].set_ylabel('x')
axes[1].set_xlabel('n')
axes[1].set_ylabel('y')
plt.show()

#最后将x和y根据正态分布函数可视化
countX, binsX, ignoredX = plt.hist(x, 30, density=True)
plt.plot(binsX, 1/(stdX * np.sqrt(2 * np.pi)) *
               np.exp( - (binsX - meanX)**2 / (2 * stdX**2) ),
         linewidth=2, color='r')
plt.xlabel('x')
plt.show()

countY, binsY, ignoredY = plt.hist(y, 30, density=True)
plt.plot(binsY, 1/(stdY * np.sqrt(2 * np.pi)) *
               np.exp( - (binsY - meanY)**2 / (2 * stdY**2) ),
         linewidth=2, color='b')
plt.xlabel('y')
plt.show()

由此可以看出,在实际计算当中,如果x是一个符合正态分布概率密度曲线的变量,其参与的计算结果也会符合相应的概率密度函数.而这种自变量的变化也会反映在因变量的计算结果当中.我们文中举的例子是一个很简单的函数关系,而实际科研当中,函数关系往往更加复杂,因此在计算当中切记要考虑数据的范围变化.
最后希望各位可以举一反三, 更好的将统计学知识应用的到科研工作中!

欢迎围观公众号 迈特派 分享更多关于科研工作中的感想和方法! 转载请注明出处