首页 杂项

菜狗完成数理统计第一次大作业的记录，记录了在完成的过程中对已学习的数理统计知识进行的回顾和总结，其实就是怕忘了 QwQ，有问题的话及时喷我，谢谢各位大佬 Orz

# 大作业内容

用 Monte Carlo 方法进行统计量分布和分位数计算

# 题目

已知统计量 T(z1,z2,...,zn)=1n∑i=1n(zi−zˉ)2zˉ−z(1)T \mathrm{(z_1,z_2,...,z_n)}=\dfrac{\sqrt[]{ {1\over n } \sum\nolimits_{i=1}^{n}(z_i-\bar{z})^2 } }{\bar z -{z_{(1)} } }，其中z1,z2,...,zn\mathrm{z_1,z_2,...,z_n} 独立同分布，均服从
（a）标准正态分布
（b）参数为 1 的指数分布

请在 (a)，(b) 两种条件下分别计算：

  1. 画出概率分布 P{T(z1,z2,...,zn)<t}P\{T\mathrm{(z_1,z_2,...,z_n)} < t\} 的近似分布曲线
  2. 对给定的 α=0.01,0.05,0.1\alpha=0.01, 0.05, 0.1，计算满足 P{T(z1,z2,...,zn)<tα}=αP\{T\mathrm{(z_1,z_2,...,z_n)} < t_\alpha\}=\alpha 的分位数 tαt_\alpha。

# 目的

提升统计软件的使用能力（本题目可以使用任意软件），熟练使用 Monte Carlo 方法进行概率计算，并能把其用于科学研究。

# 基本要求

  1. 用 Monte Carlo 方法给出计算方法过程，输出计算结果。
  2. 讨论 Monte Carlo 方法计算结果随着样本容量变化时计算误差的变化，并进行总结。
  3. 作业书写图文并茂，具有较强的可读性与可视性。

注：除完成基本要求外，可以增加自己感兴趣的研究内容。

# 算法引入

考虑这么一种情况：我们在一张纸上画了一个圆圈，此时我们想求圆圈的面积。如果圆圈时正圆形的话，面试非常好求，但如果不是一个正圆，该怎么求面积呢？

Mentor Carlo问题引入

此时我们可以考虑这么一种情况：抓一把绿豆，随机的撒在整个纸上，然后去计算 NcircleNpaper×Spaper\frac{N_{circle} }{N_{paper} } \times S_{paper} 即可计算出圆的面积。显然的是，绿豆的个数越多，算出的圆的面积就会越接近实际答案，这其实就是 蒙特卡洛（Mento Carlo） 算法的基本思想。

从上述可以看出，Monte-Carlo 算法区别于确定性算法，它的解不一定是准确或正确的，其准确性或正确性依赖于概率和统计，但在某些问题上，当重复实验次数足够大时，可以从很大概率上（这个概率是可以在数学上证明的，但依赖于具体问题）确保解的准确或正确性，所以，我们可以根据具体的概率分析，设定实验的次数，从而将误差或错误率降到一个可容忍的程度。

下面给出严格的证明：

证明

为不失一般性，可以设纸的面积为 SS，不规则图形面积为 S′S' ，设事件 A：随机投掷一粒绿豆，绿豆落在圆内。

显然，投掷落点符合二维均匀分布，所以 p(A)=S′Sp(A)=\frac{S'}{S}。

设 kk 是 nn 次投掷中，都在不规则圆形的次数，ε>0\varepsilon>0 为任意正数，根据伯努利大数定律：

lim⁡n→∞P{∣kn−p(A)∣<ε}=lim⁡n→∞P{∣kn−S′S<ε∣}=1\lim_{n \to \infty} P\left \{ \left | \frac{k}{n} - p(A) \right |< \varepsilon \right \} = \lim_{n \to \infty} P \left \{ \left | \frac{k}{n} - \frac{S'}{S} < \varepsilon \right | \right \} = 1

这就证明了，当 n→∞n \rightarrow \infty 时，频率 knk\over n 依概率收敛于 S′SS'\over S, 由于纸的面积是很好算的，也就可以轻松算出不规则圆的面积了。

上述证明从数学上说明用频率估计不规则图形面积的合理性，进一步可以给出误差分析，从而选择合适的实验次数 nn，来将误差控制在可以容忍的范围内，Monte-Carlo 算法虽然不能保证解一定是准确和正确，但并不是 “撞大运”，其正确性和准确性依赖概率论，有严格的数学基础，也因此 Mento Carlo 算法也称统计模拟法、统计实验法。

# 解题流程

# 画出概率分布近似分布曲线

由于是使用 Monte Carlo 算法来近似的求解问题，那么我们首先需要生成足够多的满足题目要求的随机数据，这里以 “标准正态分布” 数据为例进行说明：

首先我们要生成足够多的数据去 “撒绿豆”，先随机生成一组样本容量为 2000 的标准正态分布样本数据，并画出他的概率密度曲线：

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
num_arr = np.random.normal(0,1,2000)
sns.histplot(num_arr, kde=True)
plt.show()
标准正态分布图
根据 T(z1,z2,...,zn)T \mathrm{(z_1,z_2,...,z_n)} 的公式算出一个 tt，并将其加入到：
t_arr = []
z_a = np.mean(num_arr)
z_1 = np.min(num_arr)
n = len(num_arr)
tmp = 0.0
for i in range(n):
  tmp += (num_arr[i]-z_a)**2
t = math.sqrt(tmp/n) / (z_a - z_1)
t_arr.append(t)
此时，只需要重复这个算法，算出多个 tt 值去估计概率分布即可。但是还有一个需要注意的点，由于 seaborn 库函数的原因，如果想画出概率密度的话，必须将生成的数组格式，使用 pandas 库转换成 DataFrame 格式才可以正确处理。
最终代码如下：
import math
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import pandas as pd
sns.set_style("whitegrid")
t_arr = []
loop = 1000
for _ in range(loop):
    num_arr = np.random.normal(0,1,1000)
    z_a = np.mean(num_arr) # 样本平均值
    z_1 = np.min(num_arr) # 样本最小值
    n = len(num_arr) # 样本容量
    tmp = 0.0
    for i in range(n):
        tmp += (num_arr[i]-z_a)**2 # 中心二阶距
    t = math.sqrt(tmp/n) / (z_a - z_1) # t 值
    t_arr.append(t)
t_arr.sort(reverse=True) # 对 t 进行排序，方便计算概率分布
t_arr_final = []
for k in range(loop): # 格式化列表， 为生成 DataFrame 做准备
    if ( k < (loop - 1) and t_arr[k] > t_arr[k+1]):
        t_arr_final.append([t_arr[k], (loop-k)/loop])
    elif ( k == loop - 1 ):
        t_arr_final.append([t_arr[k], 0])
    elif (t_arr[k] == t_arr[k+1]):
        continue
data = pd.DataFrame(data=t_arr_final, columns=['value', 'probability'])
sns.scatterplot(x='value', y='probability', data=data) # 绘制散点图
sns.lineplot(data=data, x='value', y='probability',estimator='max', color='red') # 绘制分布曲线
plt.show() # 展示
函数概率分布
当 z1,z2,...,znz_1, z_2, ..., z_n 服从参数为 1 的指数分布时，整体代码不变，只需要将随机值生成时用到的函数改变一下，使生成的随机值满足指数分布即可，代码只需要变动一行：
import math
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import pandas as pd
sns.set_style("whitegrid")
t_arr = []
loop = 1000
for _ in range(loop):
    # num_arr = np.random.normal(0,1,1000)
    num_arr = np.random.exponential(1, 1000) # 使随机值服从指数分布
    z_a = np.mean(num_arr) # 样本平均值
    z_1 = np.min(num_arr) # 样本最小值
    n = len(num_arr) # 样本容量
    tmp = 0.0
    for i in range(n):
        tmp += (num_arr[i]-z_a)**2 # 中心二阶距
    t = math.sqrt(tmp/n) / (z_a - z_1) # t 值
    t_arr.append(t)
t_arr.sort(reverse=True) # 对 t 进行排序，方便计算概率分布
t_arr_final = []
for k in range(loop): # 格式化列表， 为生成 DataFrame 做准备
    if ( k < (loop - 1) and t_arr[k] > t_arr[k+1]):
        t_arr_final.append([t_arr[k], (loop-k)/loop])
    elif ( k == loop - 1 ):
        t_arr_final.append([t_arr[k], 0])
    elif (t_arr[k] == t_arr[k+1]):
        continue
data = pd.DataFrame(data=t_arr_final, columns=['value', 'probability'])
sns.scatterplot(x='value', y='probability', data=data) # 绘制散点图
sns.lineplot(data=data, x='value', y='probability',estimator='max', color='red') # 绘制分布曲线
plt.show() # 展示
其概率分布如下：
指数分布的概率分布
# 求分位数
根据分位数的定义：
分位数定义
设随机变量 XX 的分布函数为 F(x)F(x)，对任意给定的实数 p(0<p<1)p(0<p<1)，若存在 xpx_p 使得
P(X⩽xp)=F(xp)=pP(X \leqslant x_p ) = F(x_p) = p
则称 xpx_p 为此概率分布的 pp 分位数。
所以其实本质上就是在分布函数上求某一个 tt，使得 F(t)=tαF(t)=t_\alpha，而这一步我们在完成列表向 DataFrame 转化的时候就已经完成了，我们只需要把他找出来就可以了。
代码如下：
# 计算分位数
a_1 =  data[data['probability'].isin([0.1])]['value'].values[0]
a_2 =  data[data['probability'].isin([0.05])]['value'].values[0]
a_3 =  data[data['probability'].isin([0.01])]['value'].values[0]
print('a = 0.1 的分位数为：' + str(a_1))
print('a = 0.05 的分位数为：' + str(a_2))
print('a = 0.01 的分位数为：' + str(a_3))
样本为标准正态分布时的分位数
同样的代码，也可以求出样本服从指数分布时的分位数：
样本为指数分布时的分位数
# 思考题
Monte Carlo 方法计算结果随着样本容量变化时计算误差的变化
我们分别给出模拟次数为 100、200、1000、10000 次时的概率分布函数曲线：
不同次数的蒙特卡洛对比图
我们可以明显看到随着模拟次数的增多，概率分布曲线越来越光滑，意味着其准确度和正确性都有较大提升，这一点也是非常直观的，我们可以考虑更极端一点的，如果只有 10 次模拟次数呢？
10模拟次数情况下的概率分布曲线
可以看到，每次的近似概率分布曲线差别非常大，这是符合我们的预期的，当模拟次数较少时，得到的结果的准确度和正确性都有较大误差，无法真实反应统计情况。
# 总结
在深入了解 Monte Carlo 算法后，我有种深深的感触：Monte Carlo 算法简直就是为计算机而生的。计算机解决问题时具有天然的计算优势，对于人类来说，解决问题的过程是拆解问题的过程。通过将复杂问题拆解成一个一个的简单问题，将不能 “计算” 的问题转化成能 “计算” 的问题，也因此，在这种过程中，人类会尽量寻找简单的模型去描述事物，哪怕与实际有些 “失真”，也是可以接受的。计算机则恰恰相反，计算机强在计算能力，强在运算速度，这时，如果还用人类解决问题的方法来让计算机解决问题，那其实就无法发挥出计算机的优势。
Monte Carlo 算法就从根本上改变了这个问题，就以求面积问题来说，我们无需输入复杂的公式，构造复杂的循环、判断逻辑来计算面积的确切值，只需要使用足够多次的模拟值，就可以以极大把握接受算出的结果，就从这次大作业来看，模拟 1 万次操作也只需 2 秒左右的时间，可谓是找到了计算机解决问题的正确逻辑。
在学习 “机器学习” 的内容时，了解到机器学习的底层逻辑是基于概率的，那时我还非常不理解为什么依赖是 “概率”，现在想来，统计与概率才是计算机的底层 “基因”，这才是计算机最擅长的东西，感谢这次大作业，感谢 Monte Carlo 算法。
# 参考
  • https://www.cnblogs.com/xumaojun/p/8541637.html
  • 孙海燕，周梦等 著，应用数理统计。北京航空航天大学出版社，2008.8
# 附录
大作业整体代码如下：
import math
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import pandas as pd
sns.set_style("whitegrid")
t_arr = []
loop = 10
for _ in range(loop):
    num_arr = np.random.normal(0,1,1000)
    # num_arr = np.random.exponential (1, 1000) # 使随机值服从指数分布
    z_a = np.mean(num_arr) # 样本平均值
    z_1 = np.min(num_arr) # 样本最小值
    n = len(num_arr) # 样本容量
    tmp = 0.0
    for i in range(n):
        tmp += (num_arr[i]-z_a)**2 # 中心二阶距
    t = math.sqrt(tmp/n) / (z_a - z_1) # t 值
    t_arr.append(t)
t_arr.sort(reverse=True) # 对 t 进行排序，方便计算概率分布
t_arr_final = []
for k in range(loop): # 格式化列表， 为生成 DataFrame 做准备
    if ( k < (loop - 1) and t_arr[k] > t_arr[k+1]):
        t_arr_final.append([t_arr[k], (loop-k)/loop])
    elif ( k == loop - 1 ):
        t_arr_final.append([t_arr[k], 0])
    elif (t_arr[k] == t_arr[k+1]):
        continue
data = pd.DataFrame(data=t_arr_final, columns=['value', 'probability'])
sns.scatterplot(x='value', y='probability', data=data) # 绘制散点图
sns.lineplot(data=data, x='value', y='probability',estimator='max', color='red') # 绘制分布曲线
plt.show() # 展示
# 计算分位数
a_1 =  data[data['probability'].isin([0.1])]['value'].values[0]
a_2 =  data[data['probability'].isin([0.05])]['value'].values[0]
a_3 =  data[data['probability'].isin([0.01])]['value'].values[0]
print('a = 0.1 的分位数为：' + str(a_1))
print('a = 0.05 的分位数为：' + str(a_2))
print('a = 0.01 的分位数为：' + str(a_3))

  
 更新于   阅读次数  次
请我喝[茶]~(￣▽￣)~*
Gality 微信支付
微信支付
Gality 支付宝
支付宝