数学理论—— 蒙特卡洛近似

1. 圆周率估算

1.1 理论

• 边长为2的正方形的横坐标范围为[-1,1]，纵坐标为[-1,1]。
• 数据点（x，y）的横坐标从[-1,1]中均匀抽样得到，纵坐标从[-1,1]中均匀抽样得到，则数据点落在圆内的概率为：$$p=\frac{A_2}{A_1}=\frac{\pi }{4}$$
• 计算误差为：$$o(\frac{1}{\sqrt{n}})$$
  
  则计算圆周率的流程为：

1. 设定一个大数n，计数器m。
2. for i = 1 to n：$$\begin{gathered} x\leftarrow [-1,1] \\ y\leftarrow [-1,1] \\ m\leftarrow m+1(当x^2+y^2\le 1时) \end{gathered}$$
3. $$\pi \leftarrow \frac{4m}{n}$$

1.2 代码实现

import random
n = 10000000
m = 0
for i in range(n):
    x = random.uniform(-1,1)
    y = random.uniform(-1,1)
    if x**2 + y**2 <= 1:
        m += 1
print(4*m/n)    

3.1415656

2. 近似积分

2.1 一元积分（univariate）

• 计算如下积分:$$I={\int }_a^{b}f(x)d_x$$

1. 从[a,b]中进行n次抽样，得到:$$x_1,x_2,...,x_n$$
2. 计算：$$Q_n=(b-a)\cdot \frac{1}{n}\cdot \sum _{i=1}^{n}f(x_i)$$
3. 可用Qn去近似积分I，因为大数定律可以保证其误差反比于：$$\sqrt{n}$$
   代码为：

import random
def f(x):
    return x**2
a = -1
b = 1
n = 100000000
sum_num = 0
for i in range(n):
    x = random.uniform(a,b)
    sum_num += f(x)
Q = (b-a)* sum_num / n 
print(Q)

结果为：

0.6666073503023225

2.2 多元积分（multivariate）

求取积分：$$I={\int }_{\Omega }f(X)d_X$$

1. 从$$\Omega$$中进行n次抽样，得到:$$X_1,X_2,...,X_n$$
2. 计算：$$V={\int }_{\Omega }{d}_X$$
3. 计算：$$Q_n=V\cdot \frac{1}{n}\cdot \sum _{i=1}^{n}f(X_i)$$
4. 可用Qn去近似积分I，因为大数定律可以保证其误差反比于：$$\sqrt{n}$$
   求取下式积分：$$\begin{gathered} f(x,y)=\{\begin{array}{rl}1& ,x^2+y^2\le 1\\ 0& ,otherwise\end{array} \\ \Omega =[-1,1]\cdot [-1,1] \\ I={\int }_{\Omega }f(x,y)d_x{d}_y \end{gathered}$$
   代码为：

def f(x,y):
    if x**2 + y**2 <= 1:
        return 1
    return 0
a1 = -1
a2 = -1
b1 = 1
b2 = 1
n = 10000000
sum_num = 0
for i in range(n):
    x = random.uniform(a1,b1)
    y = random.uniform(a2,b2)
    sum_num += f(x,y)
V = 2*2 
Q = V*sum_num / n 
print(Q)

3.1415416

3 期望估计（Estimate of Expection）

• X为d维随机变量
• P（X）为X的概率密度函数：$$P(X)\leftarrow probabilitydensityfunction$$
• $${\int }_{R}P(X)d_X=1$$
  计算f(X)的期望：

1. 依据P(X)进行随机抽样，得到：$$X_1,X_2,...,X_n$$
2. 计算：$$Q_n=\frac{1}{n}\cdot \sum _{i=1}^{n}f(X_i)$$
3. $$Q_n\approx E_{X\sim P}[f(X)]$$

代码为：

n = 10000000
sum_num = 0
for i in range(n):
    x = random.normalvariate(2,1)
    sum_num += x
print(sum_num/n)

2.000336379227587

本文部分内容为参考此视频。

by CyrusMay 2022 04 04

生命是华丽错觉
时间是贼偷走一切
————五月天（如烟）————