概统大作业-draft

摘要: 大数定律在计算机科学中的一个重要应用是随机模拟，也称蒙特卡洛算法。蒙特卡罗积分（Monte Carlo integration）是一种使用随机数进行数值积分的技术。它是一种特殊的蒙特卡罗方法，可实现对定积分进行数值计算。本文将对蒙特卡洛算法的计算原理进行阐述，探讨在函数积分估计中如何利用蒙特卡洛方法生成符合均匀分布的随机数，并通过数值模拟和理论分析验证该方法的有效性。

关键词: 大数定律；蒙特卡洛算法；计算机模拟；定积分；随机数。

1 引言

1.1问题背景

​ 在处理诸如 $I = \int_0^1 \sqrt{1-\sin{x^2}} \mathrm dx$ 的积分时，当目标函数无法不属于初等函数，无法写成出等形式，或者积分区间越来越复杂，无法直接计算出其积分值时，我们需要考量新的积分方法。这时，我们可以引入蒙特卡洛算法。

​ 蒙特卡洛方法是一种基于概率的数值计算方法，被广泛应用于估计复杂函数的积分、优化问题和随机模拟等领域，。其核心思想是通过大量随机样本来逼近期望值，进而估计积分。为了确保结果的准确性和可靠性，生成高质量的随机数是关键步骤之一。通常，我们需要生成符合均匀分布的随机数，并在此基础上进行积分估计。

​ 本文旨在探讨蒙特卡洛方法在估计函数积分时的具体实现和理论依据，并通过数值模拟和理论分析相结合的方法验证其有效性。

1.2 研究意义

​ 蒙特卡洛方法在现代科学与工程中的应用越来越广泛，其独特的随机采样策略使其在处理复杂问题、特别是高维积分和复杂系统模拟方面具有显著优势。传统数值积分方法在处理高维积分和复杂函数时，往往面临计算量巨大、误差累积等问题。蒙特卡洛方法通过随机采样，可以有效降低计算复杂度，提高计算效率。

​ 蒙特卡洛方法作为一种基于概率统计的数值方法，因其在高维积分和复杂系统模拟中的显著优势，广泛应用于物理学、金融工程、生物统计等领域。

2 蒙特卡洛积分法原理

2.1 引例

​ 首先探讨上述提及积分 $I = \int_0^1 \sqrt{1-\sin{x^2}} \mathrm dx$ 。通过简单分析可知，该积分无解析解。为了得到该积分的估计值，引入 $ Y =\sqrt{1-\sin{X^2}}$, 并使$X$服从均$[0,1]$的匀分布。因此，$X$的概率密度函数函数为：

$$f(x)= \begin{cases}1,&0\leqslant x\leqslant 1 \\0,&others \end{cases}$$

​ $Y$的均值为：

$$E(Y) = \int_0^1\sqrt{1-\sin{x^2}} f(x) \mathrm dx =I$$

​ 由此，我们可以得到原积分的积分值与Y的均值相等。为了得到$E(Y)$，我们考虑矩估计的方法，即根据Khinchine大数定律，在样本数量足够多时，我们可以将样本的算数平均值作为$Y$的数学期望：

$$\lim_{x\rightarrow+\infty}P\left(\left|\frac{1}{n}\sum_{k = 1}^n Y_k-E(Y)\right|\geqslant\varepsilon\right) = 0$$

​ 因此，只需要若干服从均匀分布的样本数据，计算其样本均值可以得到总体的矩估计值，进而得到积分值的估计值。通过计算机随机生成随机数字进行模拟，可以得到上述积分的估计值为0.8121，与实际值（0.8116）十分接近。随机数的抽取原理将在随后进行详细的探讨。

2.2 一般性方法

​ 上述引例讨论了当随机变量$X$服从均匀分布时，原积分的积分值与$Y$的均值相等。接下来，我们对于蒙特卡洛方法进行一般性的分析。

​ 对于积分 $I = \int_a^bf(x)\mathrm d x$，随机变量$X$。令$X$概率密度为$g(x), x\in[a,b]$，设$X_1,X_2,…,X_n $是一个独立同分布的随机变量序列。则可以证明得到：

$$\hat{I}_n = \frac{1}{n}\sum_{i = 1}^n \frac{f(X_i)}{g(X_i)}$$

证明如下：

$$\begin{align}I &= \int_a^bf(x)\mathrm dx \\&= \int_a^b\frac{f(x)}{g(x)}g(x)\mathrm dx \\&= E\left(\frac{f(X)}{g(X)}\right)\end{align}$$

​ 每次采样相互独立，则$\frac{f(X_i)}{g(X_i)}（i=1,2,…,n）$服从独立同分布且数学期望存在。根据Khinchine大数定律，该独立同分布的序列服从大数定律。即：

$$\forall \varepsilon >0,\lim_{x\rightarrow+\infty}P\left(\left|\frac{1}{n}\sum_{i = 1}^n\frac{f(X_i)}{g(X_i)}-E\left(\frac{f(X)}{g(X)}\right)\right|\geqslant\varepsilon\right) = 0$$

​ 令$E\left(\frac{f(X)}{g(X)}\right) = \mu$，由矩估计原理可知，$\hat{\mu} =\frac{1}{n}\sum_{i = 1}^n\frac{f(X_i)}{g(X_i)} $。即：

$$\hat{I}_n = \frac{1}{n}\sum_{i = 1}^n \frac{f(X_i)}{g(X_i)}$$

$Q.E.D$

​ 由此，我们得到了蒙特卡洛积分估计的一般方法：

1. 通过任意概率密度函数给$g(x)$产生 𝑛 个独立同分布的随机数 $X_1,X_2,…,X_n$.
2. 用这些随机数采样 $h(x)= \frac{f(x)}{g(x)}$，得到得到独立同分布的样本$\frac{f(X_i)}{g(X_i)}（i=1,2,…,n）$
3. 用$\frac{1}{n}\sum_{i = 1}^n\frac{f(X_i)}{g(X_i)} $的观测值作为$E(h(X))$ 的估计值，也即$I_n$​的估计值。

​ 蒙特卡洛方法的概率密度函数是任选的，最简单的，我们可以用均匀分布$ 𝑋∼𝑈(𝑎,𝑏) $进行估计。

2.3 随机数抽取原理

​ 对于一般的均匀分布，可以使用计算机中自带函数实现随机数生成，在现代的编程语言中，都提供了伪随机生成的库函数，比如 C++ 的 <random> 库函数。这种算法得到的“随机数列”存在一个周期，其大小为 T = 0x7fff，并且在一个周期内，每个数只出现一次，频率为 1/T 。也就是说这个随机数生成算法得到的伪随机数 𝑋 服从： $𝑋∼𝑈(0,32767)$ 。除此之外，也可以使用梅森旋转算法，对传统 PRNG 算法进行改进，高效生成高质量的服从均匀分布的伪随机数。

​ 通过 PRNG，我们能生成服从均匀分布的伪随机数。但在蒙特卡洛积分中，均匀分布并不是最合适的分布，我们需要生成服从指定分布的伪随机数。如果要生成符合条件的随机数，也有多种不同的方法可以实现。在此，笔者将引入两种常见的生成方法。

2.3.1 逆变换法

​ 使用逆变换法，我们将均匀分布的随机数映射到目标分布。具体原理如下：

$$\begin{align}设Y &=F_X(X)\\F_Y(y) &= P(Y\leqslant y) \\&= P(F(X)\leqslant y) \\&= P(X\leqslant F^{-1}(y)) \\&=F(F^{-1}(y)) \\&= y\end{align}$$

​ 因此，$Y$服从均匀分布 $Y\sim U(0,1)$.

​ 由此，我们可以使用计算机模拟出一系列服从均匀分布的随机数据，再通过逆变换得到一系列我们需要的变量的随机数据。 需要注意的是，此方法只有在分布函数可以求出反函数时有效。

示例如下：

​ 指数分布$X\sim E(1)$，其概率密度函数为：

$$f(x) =\begin{cases}e^{-x},&x>0\\0,&others \end{cases}$$

​ 分布函数为：

$$F(x) =\begin{cases}1-e^{-x},&x>0\\0,&others \end{cases}$$

​ 通过逆变换，我们求出$F(x)$​​反函数，可以得到：

$$x =F^{-1}(u) = -\ln(1-u),u\sim U(0,1)$$

​ 为验证随机数生成效果，使用python中的numpy库以及matplotlib库模拟随机抽取10000次的结果（代码见附页），并绘制频率分布直方图，将结果与理论概率密度函数作对比，可以发现得到的结果与理论值拟合情况非常好，证明模拟抽取有效。

2.3.2 取舍法

​ 设我们需要生成的随机数概率密度为$f(x)$，其中$x \in [x_{min},x_{max}],f(x) \in (0,y_{max}]$均匀分布的样本为$u$.则对于每一个随机数x，我们再通过另一个均匀分布$U(0,y_{max})$得到另一个随机数y。如果$y\leqslant f(x)$​，则x为满足条件的一个随机数。

示例如下：

​ 对于标准正态分布$X\sim N(0,1)$，其概率密度函数为：

$$f(x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} \\y_{max} = \frac{1}{\sqrt{2\pi}}$$

​ 为验证该方法是否有效，我们使用python的matplotlib和numpy库生成10000个随机样本，将样本进行取舍处理后，绘制出处理后的频率分布直方图（代码见附页），并将结果与理论概率密度函数作对比，可以发现得到的结果与理论值拟合情况较好，证明模拟抽取有效。

3 估计效果分析

3.1 理论分析

3.1.2 无偏性分析

​ 由上述推导可知，

$$\hat{I}_n = \overline{h(X)}=\frac{1}{n}\sum_{i = 1}^n \frac{f(X_i)}{g(X_i)}$$

​ 现讨论其无偏性：

$$\begin{align} \hat{I}_n &= E(\overline{h(X)} )\\&= \frac{1}{n}E\left(\sum_{i = 1}^n \frac{f(X_i)}{g(X_i)}\right) \\&= E\left( \frac{f(X)}{g(X)}\right) \\&= \int_a^b \frac{f(x)}{g(x)} \cdot g(x) \mathrm d x \\&= I_n \end{align}$$

​ 因此，样本均值$\overline{h(X)} $是$I_n$的无偏估计量。

3.1.2 有效性分析

​ 讨论其有效性，只需要计算估计量的方差：

$$D(\overline{h(X)}) = D\left(\frac{1}{n}\sum_{i = 1}^n \frac{f(X_i)}{g(X_i)}\right) =\frac{1}{n}D\left( \frac{f(X)}{g(X)}\right)$$

​ 根据方差计算公式$D(X) = E(X^2)-(E(X))^2$,

$$E\left(\left( \frac{f(X)}{g(X)}\right)^2\right) = \int_a^b \frac{f^2(x)}{g^2(x)} \cdot g(x) \mathrm dx = \int_a^b \frac{f^2(x)}{g(x)}\mathrm dx$$

​ 由于是无偏估计量，$E\left( \frac{f(X)}{g(X)}\right) = I_n$

$$D(\overline{h(X)}) = \frac{1}{n} \left(\int_a^b \frac{f^2(x)}{g(x)}\mathrm dx - I_n^2 \right)$$

​ 因此，可以得到结论：

$$D(\overline{h(X)}) \propto \frac{1}{n}$$

​ 这个式子说明，该估计量的有效性与样本数量正相关，与积分值$\int_a^b \frac{f^2(x)}{g(x)}\mathrm dx$负相关。因此，可以通过增大样本数量来增大蒙塔卡洛估计方法的精度，比如，在计算机中可以尽可能使用更多的随机数生成样本数据，从而得到更加精确的积分估计值。

3.1.3 一致性分析

​ 有上述推导可以知道，估计量$\overline{h(X)}$无偏，且，估计量的方差与样本容量成反比关系。且有：

$$\lim_{n\rightarrow+\infty}D(\overline{h(X)}) = 0$$

​ 因此，由已知结论可以知道，该估计量也是$I_n$的一致估计量。

3.2 数值模拟

​ 有上述推导可以知道，对于同一个积分，积分估计的精度与两个方面因素有关，即样本数量与积分值$\int_a^b \frac{f^2(x)}{g(x)}\mathrm dx$。接下来，本文将通过对于引例中的积分$I = \int_0^1 \sqrt{1-\sin{x^2}} \mathrm dx$，选取服从不同的概率密度的$X$，以及不同的样本数量，探讨蒙特卡洛积分方法的估计效果。

3.2.1 X分布的影响

​ 运用逆变换法，生成一系列随机样本，分别服从均匀分布、指数分布、正态分布，在[0,1]运用蒙特卡洛积分法计算积分估计值。为控制变量，样本数量设为10000.

​ 需要注意的是，在寻找满足条件的分布时，应该使$\int_0^1g(x)\mathrm d x = 1$

（1）均匀分布：

​ 随机变量满足$X$服从均匀分布$ g(x)= \begin{cases}1,&0\leqslant x\leqslant 1 \0,&others\end{cases} $

​ 使用python在[0,1]生成10000个随机数据，服从均匀分布。绘制频率直方图，如下：

​ 计算样本平均值$\frac{1}{n}\sum_{i = 1}^n \frac{f(X_i)}{g(X_i)} =\frac{1}{n}\sum_{i = 1}^n \frac{\sqrt{1-\sin{X_i^2}}}{1} $得到积分的估计值。采样五次后的数据为：

​ 0.81165；0.80939；0.81111；0.81115；0.81096

​ 样本平均值为 0.81085

（2）$g(x) =\begin{cases}2x,&0\leqslant x\leqslant 1\\0,&others \end{cases}$

​ 使用python在[0,1]生成10000个随机数据，逆变换处理后服从2x的分布。绘制频率直方图，如下：

​ 计算样本平均值$\frac{1}{n}\sum_{i = 1}^n \frac{f(X_i)}{g(X_i)} =\frac{1}{n}\sum_{i = 1}^n \frac{\sqrt{1-\sin{X_i^2}}}{2x} $得到积分的估计值。采样五次后的数据为：

​ 0.81847；0.83682；0.80306；0.79894；0.83263

​ 样本平均值为 0.81798

（3）$g(x) =\begin{cases}\frac{\pi}{2}\sin\frac{\pi}{2}x,&0\leqslant x\leqslant 1\\0,&others \end{cases}$

​ 使用python在[0,1]生成10000个随机数据，逆变换处理后,，可以得到$X = \frac{2}{\pi}\arccos(1-u),u\sim U(0,1)$.服从上述的分布。绘制频率直方图，如下：

​ 计算样本平均值$\frac{1}{n}\sum_{i = 1}^n \frac{f(X_i)}{g(X_i)} =\frac{1}{n}\sum_{i = 1}^n \frac{\sqrt{1-\sin{X_i^2}}}{\frac{\pi}{2}\sin\frac{\pi}{2}x} $​得到积分的估计值。采样五次后的数据为：

​ 0.82288；0.83586；0.83342；0.79847；0.84108

​ 样本平均值为 0.82634

（4）总结

	1	2	3	4	5	average
uniform	0.81165	0.80939	0.81111	0.81115	0.81096	0.81085
2x	0.81847	0.83682	0.80306	0.79894	0.83263	0.81798
cos	0.82288	0.83586	0.83342	0.79847	0.84108	0.82634

​ 根据相关资料，所求积分的精确值约为0.81157，因此可以发现均匀模拟最有效，精度最高。

​ 根据已有的推导可知，估计量的有效性与积分$\int_a^b \frac{f^2(x)}{g(x)}\mathrm dx$负相关。分别计算上述三个分布的对应积分值

$$\int_0^1 \frac{1-\sin^2(x)}{1}\mathrm dx = 1-\dfrac{\sqrt{\pi}\,\operatorname{S}\left(\dfrac{\sqrt{2}}{\sqrt{\pi}}\right)}{\sqrt{2}}\approx 0.6897317 \\\int{\dfrac{1-\sin\left({x}^{2}\right)}{2\,x}}{\;\mathrm{d}x}\rightarrow +\infty \\ \int_0^1{\dfrac{2\,\left(1-\sin\left({x}^{2}\right)\right)}{\pi\,\sin\left(\frac{\pi\,x}{2}\right)}}{\;\mathrm{d}x} \rightarrow +\infty$$

​ 因此，上述三种方法中，均匀分布最有效。与数值模拟的结果一致

3.2.2 样本容量的影响

​ 为了探究样本容量对于模拟计算精确的的影响，本文选取均匀分布计算上述引例的积分值，探究积分估计值与样本容量的影响。

​ 在样本数量为10，100，1000，10000，100000分别采样，将其结果与精确作差，得到误差曲线，如图所示

​ 数值模拟结果证明误差大小与样本容量成负相关关系，即样本容量越大，蒙特卡洛积分算法估计值越准确。这与之前的理论推导结果一致。同时，由于估计量方差与样本容量的一次方成反比，蒙特卡洛方法的收敛速度相对较慢，因此可能需要大量的样本数量才能达到较高的精度。

​ 由此回归引例，可以在尽可能增大样本容量的条件下，我们给出该积分的高精度解：0.811573640825

4 方差缩减

​ **方差缩减（variance reduction）**是在不增加蒙特卡洛试验次数的前提下，减小估计的误差。被积函数f(x)是确定的，概率密度函数g(x)是任选的，因此能不能找到一个最合适的概率密度函数，使得积分值$\int_a^b \frac{f^2(x)}{g(x)}\mathrm dx$最小，进而使得方差最小？我们讨论$D(\overline{h(X)}) = \frac{1}{n} \left(\int_a^b \frac{f^2(x)}{g(x)}\mathrm dx - I_n^2 \right)$，可以得出，当：

$$g(x) = \frac{|f(x)|}{\int_a^bf(x)\mathrm d x} \\D(\overline{h(X)})_{min} = 0$$

​ 但是，f(x)的积分难以计算，那么蒙特卡洛积分又退回到了定积分的计算。但这也给我们方差缩减提供了思路。因此我们引入重要性采样。

4.1 重要性采样

​ 蒙特卡洛法的核心，是根据某一概率分布来随机采样，即根据被积函数曲线规律来设计类似的概率密度来进行采样，这就符合蒙特卡洛积分的要求，也就是重要性采样的核心。

事实上这一理想的$g(x)$不可能得到，不过我们可以选择与被积函数形状相似的$g(x)$，它可使得$\frac{|f(x)|}{\int_a^bf(x)\mathrm d x}$趋近于一个常数，这能有效的减小蒙特卡洛积分的方差。

​ 接下来以计算$\int_0^{\frac{\pi}{2}}\sin x\mathrm d x$为例，对重要性采样进行说明。

​ 我们设计两种概率密度函数，一种是均匀分布，随机变量$X_1\sim U(0,\frac{\pi}{2})$；另外一种，随机变量$X_2$概率密度函数为$g(x) =\begin{cases}\frac{8x}{\pi^2},&0\leqslant x\leqslant \frac{\pi}{2},\\0,&others \end{cases}$

​ 两种概率密度与原被积函数作比较，绘图得到：

​ 由图可知，第二种分布相比于均匀分布更加接近原被积函数，因此理论推导上第二种会更加精确。接下来，使用数值模拟验证。

（1）均匀分布：

​ 随机变量满足$X_1$服从均匀分布$ g(x)= \begin{cases}\frac{2}{\pi}&0\leqslant x\leqslant \frac{\pi}{2} \0,&others\end{cases} $

​ 使用python在[0,1]生成10000个随机数据，计算所求估计量的样本平均值，重复五次，得到数据为：

1.00596；1.00274；0.99942；0.99664；0.99406

​ 其平均值为 0.999764，残差平方和约为 8.9939E-05

（2）$X_2$:

​ 随机变量满足$X_2$服从分布$g(x) =\begin{cases}\frac{8x}{\pi^2},&0\leqslant x\leqslant \frac{\pi}{2},\\0,&others \end{cases}$

​ 使用python在[0,1]生成10000个随机数据，计算所求估计量的样本平均值，重复五次，得到数据为：

1.00168；1.00096；1.00042；0.99921；1.00067

​ 其平均值为 1.000588 ，残差平方和约为 4.9934E-06

​ 由此可以看出，在相同条件下，重要性采样估计值比均匀采样的误差小。

4.2 分层采样

​ 分层采样（Stratified Sampling）把f(x)的积分区间均匀的划分成 m个不相交的子区间$\Omega_i$称为阶层（Stratum）。在每个阶层内用$n_i$分别进行蒙特卡洛积分。整个积分域内的蒙特卡洛积分就是各个阶层蒙特卡洛积分的和。分层采样通常可以采样均匀分层。但随着被积函数维数不断增加，分出来的层数呈指数级增加，为了解决这个问题，人们提出了拟蒙特卡洛方法。

​ 拟蒙特卡洛积分估计技术的核心点是积分估计时采用低差异序列（Low Discrepancy Sequence）来替换纯随机数，它的优点是积分估计的收敛速度更快，理论上拟蒙特卡洛能达到$o(n^{-1})$的收敛速度，而普通蒙特卡洛的收敛速度是$O(n^{-0.5})$​，更快的收敛速度意味着相同的采样数下，低差异序列估计的结果误差更小。

​ 但是，分层采样方法超出了目前的知识范围，本文在此不做深究。

5 近似计算方法比较

​ 其他积分近似计算的方法有很多，比如矩形法、梯形法、辛普森法等等。本文在此以计算积分$\int_0^1x^2 \mathrm d x$​​为例，分析矩形法与蒙特卡洛方法的比较。

（1）矩形法

n	10	100	1000	10000
I	0.3325	0.333325	0.33333325	0.333333333

（2）蒙特卡洛法

​ 由重要性采样，可令$g(x) =\begin{cases}1,&0\leqslant x\leqslant 1\\0,&others \end{cases}$

​ 使用python在[0,1]生成10、100、1000、10000、100000个随机数据，经过逆变换处理后，计算所求估计量的样本平均值，得到数据为：

n	10	100	1000	10000	100000
I	0.515564	0.371139	0.335690877	0.332330614	0.333225706

​ 由此可见，矩形法的收敛速度快于蒙特卡洛法。从原理上来说，微积分近似适用于能够使用解析方法得到积分的函数，通常对于简单函数和特定形状的曲线效果较好。而蒙特卡洛积分法更灵活，适用于高维、复杂的积分问题，对于无法通过解析方法求解的问题效果更为明显。同时，蒙特卡洛法可以任选服从任意概率密度的分布，灵活性更强，编程操作更容易实现。

6 结论

​ 本文从蒙特卡洛法的原理入手，得出蒙特卡洛积分法的一般性方法，即对于积分 $I = \int_a^bf(x)\mathrm d x$，随机变量$X$。令$X$概率密度为$g(x), x\in[a,b]$，设$X_1,X_2,…,X_n $是一个独立同分布的随机变量序列。则有：

$$\hat{I}_n = \frac{1}{n}\sum_{i = 1}^n \frac{f(X_i)}{g(X_i)}$$

​ 进而，对于随机变量的抽取，本文引入了两种方法，即逆变换法和取舍法。在之后的效果分析中，笔者从理论分析的数值模拟两方面证明了蒙特卡洛方法的精度取决于概率密度函数的选取，以及样本容量的大小，并引入重要性采样的方法。最终得出结论：

• X的选取一定时，$D(\overline{h(X)}) \propto \frac{1}{n}$
• 样本容量一定时，当$g(x) = \frac{|f(x)|}{\int_a^bf(x)\mathrm d x} ，D(\overline{h(X)})_{min} = 0$

​ 蒙特卡洛模拟法，是大数定理在定积分近似计算中的一个应用，是一种灵活、适用于复杂问题的定积分近似计算方法。研究生成均匀分布随机数的方法和评估其效果，不仅能提升数值计算能力，还能优化计算精度与效率，对随机模拟技术和实际应用具有重要意义。

​ 蒙特卡洛的优越性更多体现在应对多维积分时的近似计算，本文只体现了应对一维积分的计算方法。在后续研究中，需要使用到分层采样的方法，包括拟蒙特卡洛方法和N-Rooks方法等，探究多维度的蒙特卡洛应用。这些都超出了目前的知识范围，可在以后的内容中探讨。

附录

本文用到的蒙特卡洛模拟代码如下：

#逆变换
import random
from math import sin,sqrt
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import gaussian_kde
import seaborn as sns
list1 = []
sum = 0
for _ in range(100000):
    random_number = round(random.random(),4)
    list1.append(random_number)
for i in list1:
    sum+=(sqrt(1-sin(i**2)))
y1 = np.random.rand(10000) 
y2 =-np.log(1-y1)
kde = gaussian_kde(y2)
x_values = np.linspace(min(y2), max(y2), 10000)
pdf_values = kde(x_values)
plt.hist(y2, bins=100, density=True, alpha=0.5, color='b',
label='Generated Data')
plt.scatter(x_values,pdf_values,s=10)
x_values_exp = np.linspace(min(y2), max(y2),10000)
y_values_exp = np.exp(-x_values_exp)
plt.plot(x_values_exp, y_values_exp, color='red', label='$e^{-x}$')
plt.title('exponential distribution')
plt.xlabel('X')
plt.ylabel('Y')

plt.legend()
plt.show()

#数值模拟
#uniform
import numpy as np
import matplotlib.pyplot as plt
from math import pi,sin,e,sqrt

y1 = np.random.rand(10000) 
yy1 = list(y1)
sum = 0
for i in yy1:
    sum+= sqrt(1-sin(i**2))
print(sum/10000)
plt.hist(y1, bins=50, density=True, alpha=0.5, color='b',label='Generated Data')

plt.title('uniform distribution')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.show()

#2x
from math import sin,sqrt,e
import matplotlib.pyplot as plt
import numpy as np

def f(x):
    return sqrt(1-sin(x**2))
def g(x):
    return 2*x

y1 = np.random.rand(10000)  

y2 = np.sqrt(y1)              
list2 = list(y2)
sum_list = 0
for i in list2:
    sum_list+=f(i)/(g(i))
sum_list /= 10000
print(round(sum_list,5))

print(sum/10000)
plt.hist(y2, bins=100, density=True, alpha=0.5, color='b',label='Generated Data')

plt.title('2x distribution')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.show()

#sin
from math import sin,sqrt,e,pow,pi
import matplotlib.pyplot as plt
import numpy as np

def f(x):
    return sqrt(1-sin(x**2))
def g(x):
    return pi/2*sin((pi/2)*x)

y1 = np.random.rand(10000) 
y2 = (2/pi) *np.arccos(1-y1) 
list1 =list(y2)
res = 0
for i in list1:
    res+=f(i)/g(i)
res/=10000

print(round(res,5))
plt.hist(y2, bins=100, density=True, alpha=0.5, color='b',label='Generated Data')
plt.title('sin distribution')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.show()

#重要性采样
import numpy as np
import matplotlib.pyplot as plt
from math import pi,sin,e,sqrt

def f1(x):
    return sin(x)
def g1(x):
    return 2/pi

def f2(x):
    return sin(x)
def g2(x):
    return 8*x/(pi**2)

y1 = np.random.uniform(0,1,10000) 
yy1 = list(y1)
y2 = pi/2*(np.sqrt(y1))
yy2 = list(y2)
sum = 0
for i in yy2:
    sum+= f2(i)/g2(i)
print(round(sum/10000,5))

x = np.linspace(0,pi/2,10000)
y1 = x- x+2/pi
plt.plot(x, y1, color='red', label='y1')
y2 = 8*x/(pi**2)
y3 = np.sin(x)
plt.xlabel('X')
plt.ylabel('Y')
plt.plot(x, y2, color='blue', label='y2')
plt.plot(x, y3, color='g', label='sinx')
plt.title("importance compare")

plt.legend()
plt.show()

#矩形法
import math
import numpy as np

samples = 10000
lenth = 1/samples
x =np.linspace(1/samples/2,1-1/samples/2,samples)
y = x**2
res = 0
for i in y:
    res+=i*lenth
print(res)

​