Python蒙特·卡罗方法和Q-Q plot验证中心极限定理

一.蒙特卡罗方法

蒙特卡洛的基本原理简单描述是先大量模拟，然后计算一个事件发生的次数，再通过这个发生次数除以总模拟次数，得到想要的结果，精髓就是：用统计结果去计算频率，从而得到真实值的近似值。蒙特卡洛方法可以应用在很多场合，但求的是近似解，在模拟样本数越大的情况下，越接近与真实值，但样本数增加会带来计算量的大幅上升。

不理解的话请戳:
https://blog.csdn.net/crj0926/article/details/101052873

二.Q-Q plot

• Q，Quantile，分位数，亦称分位点。

  不理解分位点戳：https://blog.csdn.net/crj0926/article/details/100853490

• 用于直观验证一组数据是否来自某个分布，或者验证某两组数据是否来自同一（族）分布。
  经常用于检验数据是否来自于正态分布。

• 用QQ图去验证数据是否满足正态分布，只需看QQ图上的点是否近似地在一条直线附近。

三.中心极限定理

四.思路

• 大概就是你先选择一种分布产生数据,我这里选择的是均匀分布(期望和方差比较好算),然后去验证产生的数据是否满足:
  
  如果满足:构建的QQ图上的点近似地在一条直线附近。
  

五.代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import matplotlib.pyplot as plt
import random
import pandas as pd
import numpy as np
from numpy import nan as NA
from scipy import stats
#####随机生成数字并取平均值####
random_data = np.random.uniform(-1, 1, 20000)#均匀生成20000个整数（在-1和1之间）
samples_mean = [] #平均值的数组
for _ in range(20000):  #这个是范围（对取到的10000个整数操作）
    sample=np.random.choice(random_data,20000)#取20000个数然后算平均值
    samples_mean.append(sample.mean())
samples_mean=np.array( samples_mean) 
plt.hist(samples_mean,bins=30,color='g')#bins=直方图中产生多少个柱状图
plt.grid()
plt.show()

q = norm.cdf(3, 0, 1) # 正态分布函数值  μ=0  σ=1
norm.ppf(q)
n = 20000  #取到值的数量
mu = 0     # 这里为取到的所有值的平均值即E(X)=0
sigma = pow(1/3, 1.0/2)  # 这里为取到的数组标准差的值,(-1,1),用公式
#利用pow(a, b)函数即可。需要开a的r次方则pow(a, 1.0/r)。
#因为这里我用到的是平均分布产生数据
#平均分布期望E(X)=(a+b)/2,方差为D(X)=(b-a)^2/12,根据你取平均值的范围求这两个值

###################
##构建Q Qplot图#####
###################
y = (samples_mean-mu)*pow(20000, 1.0/2)/sigma
y.sort()
prob = (np.arange(n)+1/2)/n
q = norm.ppf(prob, 0, 1)

plt.scatter(x = q, y = y, color = 'red')
plt.plot([y[0], y[n-1]], [y[0], y[n-1]], color = 'blue', linewidth = 2)
plt.xlabel('Theoretical Quantiles')
plt.ylabel('Ordered Values' )
plt.title('Normal Q-Q plot')
plt.show()

六.结果

示例代码输出结果:

基本满足正态分布,可以证明中心极限定理。

七.随机数证明中心极限定理

这种方法我们可以宏观的证明一下中心极限定理
首先,我们通俗的理解一下中心极限定理:
比如我们有X个数,我们把这X个数随机分成m组,然后对每一组求一个平均值n,如果我们得到无穷多个n值,它应该满足正态分布.

代码:

random_data = np.random.randint(-1, 1, 10000)#随机生成10000个整数（只有1和2）
samples_mean = []
for _ in range(10000):
    sample=np.random.choice(random_data,10000)
    samples_mean.append(sample.mean())
samples_mean=np.array( samples_mean)
plt.hist(samples_mean,bins=30,color='g')
plt.grid()
plt.show()

下面是我得到的条形图,从条形图上大概可以看出越来越满足正态分布.

2W个数每次取2000个数操作：

2W个数每次取2W个数操作：

100W个数每次取100W个数操作：

1000W个数每次取1000W个数操作：