蒙特卡洛方法原理与应用统计量近似分布与分位数计算 python画图演示

本文链接：https://blog.csdn.net/weixin_45792450/article/details/133638257

本文介绍了蒙特卡洛方法的基本原理和起源，并通过实例展示了如何利用该方法进行统计量分布和分位数计算。通过Python编程实现，探讨了样本容量变化对计算误差的影响以及分位数对随机数初始值的敏感性。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

一蒙特卡洛概率方法的初步介绍与基本原理

1.1 蒙特卡洛方法的定义

蒙特卡洛方法又称统计模拟法、随机抽样技术，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。使用蒙特卡洛方法时，需要将所求解的问题同一定的概率模型相联系，然后用电子计算机实现统计模拟或抽样，以获得问题的近似解。

1.2 蒙特卡洛方法的起源

蒙特卡洛方法于20世纪40年代美国在第二次世界大战中研制原子弹的"曼哈顿计划"计划的成员S.M.乌拉姆和J.冯·诺伊曼首先提出。数学家冯·诺伊曼用驰名世界的赌城-摩纳哥的Monte Carlo来命名这种方法，为它蒙上了一层神秘色彩。蒙特卡洛方法在金融工程学，宏观经济学，生物医学，计算物理学(如粒子输运计算、量子热力学计算、空气动力学计算、核工程)等领域应用广泛。

1.3 蒙特卡洛方法的举例

我们可以使用蒙特卡洛方法巧算 $π\pi$ 的近似值，具体过程为：正方形内部有一个相切的圆，它们的面积之比是 $π4\frac{\pi}{4}$ ，如下图所示。

此时，如果我们在正方形内部产生随机的n个点（n足够大），那么圆形内部的点数量与正方形内的点数量之比，将在概率意义上近似等于圆与正方形的面积之比，假设这个比值为 $k$ ，那么便有公式： $π=4k\pi = 4k$ 。

二用蒙特卡洛进行统计量分布和分位数计算

2.1 数理统计问题的描述

2.2 正态分布与指数分布

2.3 经验分布与概率计算

2.4 解决思路与大概过程

        我们以标准正态分布的案例进行说明，相应指数分布的思路是差不多的。在这一节，我们主要叙述大概步骤，具体的答案在第四章里面给出。
        第一步，从标准正态分布 $N (0, 1)$ 中随机取样获得足够多的 $m$ 个点： $x_{1} ,x_{2},...,x_{m}$
        第二步，将 $m$ 个点随机分成 $k$ 组，每组含 $n$ 个点，并按题中公式算出： $t_{1} ,t_{2},...,t_{k}$
        第三步，对获得的 $k$ 个点 $t_{1} ,t_{2},...,t_{k}$ ，按照经验分布定义求出对应近似分布表达式
        第四步，利用上述获得的经验分布表达式，去完成我们目标统计量分位数计算任务

三蒙特卡洛方法去解决问题的实际具体过程

3.1 相应分布随机数获取

我们可以利用python代码获取正态分布和指数分布的随机数，相关代码如下：

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

# 生成n*10万个标准正态分布随机数和标准指数分布随机数
n = 20
dataX_nor = np.random.normal(0, 1, n * 10 * 10000)
dataX_exp = np.random.exponential(1, n * 10 * 10000)

然后，我们需要对获取到的100万个随机数以每组 $n$ 个形式进行分组，相关代码如下：

# 将生成的正态分布和指数分布随机数按n个一组形式分组
dataZ_nor = np.reshape(dataX_nor, (n * 10 * 10000 // n, n))
dataZ_exp = np.reshape(dataX_exp, (n * 10 * 10000 // n, n))
# 将得到的分组随机数去计算目标统计量T的若干观测值
dataT_nor = np.zeros(n * 10 * 10000 // n, np.float32)  # 创建空数组
dataT_exp = np.zeros(n * 10 * 10000 // n, np.float32)  # 创建空数组
for i in range(n * 10 * 10000 // n):
    sum_nor, sum_exp = 0, 0
    for j in range(n):
        sum_nor += (dataZ_nor[i][j] - dataZ_nor[i].mean()) ** 2
        sum_exp += (dataZ_exp[i][j] - dataZ_exp[i].mean()) ** 2
    dataT_nor[i] = math.sqrt(sum_nor / n) / (dataZ_nor[i].mean() - dataZ_nor[i].min())
    dataT_exp[i] = math.sqrt(sum_exp / n) / (dataZ_exp[i].mean() - dataZ_exp[i].min())

3.2 近似分布函数的绘画

我们可以构造一个程序函数去代表近似分布函数 $F (x)$ ，程序函数的输入为 $x$ 和数据集，输出为 $F (x)$ ，处理过程按照经验分布的定义即可，相关代码如下：

# 利用程序去表示近似分布函数，输入自变量，输出函数值
def F(x, data):
    k = len(data)
    u = 0
    for index in range(k):
        if x > data[index]:
            u += 1
    return u / k

有了近似分布函数后，我们就可以绘画其近似分布函数的图像了，相关代码如下：

# 绘制函数图像
X = np.linspace(0, 1.5, 150)  # 创建从-1到9的等间隔的100个数据点
Y_nor = np.zeros(len(X), np.float32)  # 创建空数组
Y_exp = np.zeros(len(X), np.float32)  # 创建空数组
for index in range(len(X)):
    Y_nor[index] = F(X[index], dataT_nor)
    Y_exp[index] = F(X[index], dataT_exp)
plt.plot(X, Y_nor)  # 绘制函数图像
plt.xlabel('x')  # 设置x轴标签
plt.ylabel('y')  # 设置y轴标签
plt.title('Function-Nor')  # 设置图像标题
plt.grid(True)  # 显示网格线
plt.show()  # 显示图像

有了上述的步骤和代码，我们大致得到了如下近似分布曲线，其中标有 $N or$ 的曲线计算数据来源正态分布，标有 $E x p$ 的曲线计算数据来自于指数分布，两者均满足 $n = 10$ ：

3.3 目标统计量的分位数

从近似分布曲线图像中，我们很容易大致估计所求分位数在区间 $(0, 1)$ 中,根据经验分布和顺序统计量相关知识，我们很容易得到公式 $tα=T(kα)t_{\alpha}=T_{(k\alpha )}$ ，我们可以对 $t_{1} ,t_{2},...,t_{k}$ 排序后，按照此公式取相应的值，相关代码如下：

# 定义函数计算统计量分位数
def A(alpha, data):
    k = len(data)
    data_sort = np.sort(data)  # 排序
    return data_sort[int(k * alpha)]

取 $n = 10$ ，根据程序运行结果，我们得到：

	标准正态分布	标准指数分布
$α=0.01\alpha =0.01$	0.3936533	0.56543356
$α=0.05\alpha =0.05$	0.4370795	0.65696585
$α=0.10\alpha =0.10$	0.4659587	0.7118661

3.4 样本容量变化的误差

为了观察分位数值与样本容量关系，我们分别取值 $n = 10, n = 20, ..., n = 70$ 进行观察。
对于标准正态分布，结果如下：

$α\alpha$ \ n	$n = 10$	$n = 20$	$n = 30$	$n = 40$	$n = 50$	$n = 60$	$n = 70$
0.01	0.39365	0.33851	0.31550	0.30416	0.29674	0.29148	0.28938
0.05	0.43708	0.38108	0.35766	0.34427	0.33517	0.32763	0.32317
0.10	0.46596	0.40843	0.38308	0.36831	0.35759	0.34942	0.34489

对于参数为1的标准指数分布，结果如下：

$α\alpha$ \ n	$n = 10$	$n = 20$	$n = 30$	$n = 40$	$n = 50$	$n = 60$	$n = 70$
0.01	0.56543	0.64118	0.69019	0.71926	0.74442	0.76087	0.77441
0.05	0.65697	0.72240	0.76146	0.78683	0.80469	0.81802	0.82034
0.10	0.71187	0.76891	0.80247	0.82376	0.83893	0.84743	0.85112

通过上述的运行结果，我们容易发现，随着 $n$ 值的增大，计算结果趋向稳定，并且呈现单调递增或单调递减的规律，据此可以猜测计算误差正在变小。

3.5 分位数与随机性关系

作为题外要求的探讨，我们新增了一个研究项目：我们计算的分位数结果与初始随机数有何关系？如果初始随机数选取变化了，那么结果是否大幅变化？在这一节，我们探究分位数计算结果是否对随机数初始值敏感，为此我们设定 $n=100,α=0.05n=100,\alpha = 0.05$ ，然后不断改变初始随机数来观察计算结果的变化，相关表格如下：

	随机结果1	随机结果2	随机结果3	随机结果4	随机结果5
标准正态分布	0.31474	0.30983	0.31064	0.30978	0.31054
标准指数分布	0.84381	0.84988	0.85204	0.84550	0.85045

通过上述表格，我们容易发现，虽然计算结果有所波动，但大致都落在一定合理范围，呈现明显的稳定性。

四本篇文章相关参考资料与参考文献的声明

蒙特卡洛方法的介绍与原理来自包括百度百科等互联网资料站点
用蒙特卡洛进行统计量分布和分位数计算部分参考韦卫老师课件
编程语言为python，编辑器为pycharm，环境为anaconda+Win10
文章中的数理统计题目为本人研究生作业，里面代码为本人亲自写

整体代码如下，某些代码在获取不同结果时有轻微改动：

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

# 生成n*10万个标准正态分布随机数和标准指数分布随机数
n = 100
dataX_nor = np.random.normal(0, 1, n * 10 * 10000)
dataX_exp = np.random.exponential(1, n * 10 * 10000)
# 将生成的正态分布和指数分布随机数按n个一组形式分组
dataZ_nor = np.reshape(dataX_nor, (n * 10 * 10000 // n, n))
dataZ_exp = np.reshape(dataX_exp, (n * 10 * 10000 // n, n))
# 将得到的分组随机数去计算目标统计量T的若干观测值
dataT_nor = np.zeros(n * 10 * 10000 // n, np.float32)  # 创建空数组
dataT_exp = np.zeros(n * 10 * 10000 // n, np.float32)  # 创建空数组
for i in range(n * 10 * 10000 // n):
    sum_nor, sum_exp = 0, 0
    for j in range(n):
        sum_nor += (dataZ_nor[i][j] - dataZ_nor[i].mean()) ** 2
        sum_exp += (dataZ_exp[i][j] - dataZ_exp[i].mean()) ** 2
    dataT_nor[i] = math.sqrt(sum_nor / n) / (dataZ_nor[i].mean() - dataZ_nor[i].min())
    dataT_exp[i] = math.sqrt(sum_exp / n) / (dataZ_exp[i].mean() - dataZ_exp[i].min())


# 利用程序去表示近似分布函数，输入自变量，输出函数值
def F(x, data):
    k = len(data)
    u = 0
    for index in range(k):
        if x > data[index]:
            u += 1
    return u / k


# 定义函数计算统计量分位数
def A(alpha, data):
    k = len(data)
    data_sort = np.sort(data)  # 排序
    return data_sort[int(k * alpha)]


print(A(0.05, dataT_nor))
print(A(0.05, dataT_exp))
# 绘制函数图像
X = np.linspace(0, 2, 200)  # 创建从-1到9的等间隔的100个数据点
Y_nor = np.zeros(len(X), np.float32)  # 创建空数组
Y_exp = np.zeros(len(X), np.float32)  # 创建空数组
for index in range(len(X)):
    Y_nor[index] = F(X[index], dataT_nor)
    Y_exp[index] = F(X[index], dataT_exp)
plt.plot(X, Y_exp)  # 绘制函数图像
plt.xlabel('x')  # 设置x轴标签
plt.ylabel('y')  # 设置y轴标签
plt.title('Function-Exp')  # 设置图像标题
plt.grid(True)  # 显示网格线
plt.show()  # 显示图像