文章目录
蒙特卡洛方法(Monte Carlo Simulation),蒙特卡洛方法又称“随机抽样方法”,和一般数值计算方法有本质区别的计算方法,是一种近似推断的方法,属于试验数学的分支,利用随机数进行统计试验,以求得统计特征,比如通过采样大量粒子的方法来求解期望、均值、面积、积分等问题。
本文介绍直接采样、接受拒绝采样与重要性采样三种蒙特卡洛采样方法:
- 直接采样最简单,但是需要已知累积分布函数(CDF)f(x)的形式,并且好求逆函数,得到服从分布f(x)的样本.
- 接受拒绝采样是给出一个提议分布 p ( x ) p(x) p(x)(概率密度函数,PDF),对满足原分布的粒子进行采样,得到服从p(x)(PDF)分布的样本.
- 重要性采样则是给予每个粒子不同的权重,计算f(x)在概率密度函数为p(x)下的均值: E [ f ( x ) ] = ∫ f ( x ) p ( x ) d x \mathbb{E}[f(x)]=\int f(x)p(x)dx E[f(x)]=∫f(x)p(x)dx.
一、均匀分布采样
首先我们介绍均匀分布采样,这是其他采样方法的基础。均匀分布Uniform(0,1)
的样本是相对容易生成的,主要原理是线性同余随机数生成器。 通过线性同余发生器(LCG)可以生成伪随机数,我们用确定性算法生成[0,1]之间的伪随机数序列后, 这些序列的各种统计指标和均匀分布Uniform(0,1)的理论计算结果非常接近。这样的伪随机序列就有比较好的统计性质,可以被当成真实的随机数使用。
线性同余随机数生成器如下: x n + 1 = ( a x n + c ) m o d m x_{n+1}=(ax_n+c)\mod m xn+1=(axn+c)modm式中
a
,c
,m
是数学推导出的合适的常数。这种算法产生的下一个随机数完全依赖当前的随机数,所以当随机数序列足够大的时候,随机数也会出现重复子序列的情况,这也是为什么说计算机产生的随机数是伪随机数
。
总的来说,服从均匀分布采样的粒子我们可以很容易获得,均匀分布采样也是后面几种采样的基础。后文我们都假设能够直接从计算机获得服从均匀分布采样的粒子。
二、直接采样
直接采样的方法是根据累积分布函数进行采样,得到服从某个概率分布的样本。对一个已知概率密度函数(比如均匀分布)与累积概率密度函数的概率分布,我们可以直接从**累积分布函数(CDF)**进行采样。
原理:
例
设 x x x服从累积分布函数为 F ( x ) = 1 − e − x F(x)=1-e^{-x} F(x)=1−e−x(可验证是单调不减,且积分为1的函数)的分布,则可以通过逆变换的方法对 F ( x ) F(x) F(x)直接采样,产生服从F(X)分布的样本X.
令 y = 1 − e − x → e − x = 1 − y 两边求对数可得 : x = − l n ( 1 − y ) 令\quad y=1-e^{-x} \rightarrow e^{-x}=1-y \\ 两边求对数可得:x=-ln(1-y) 令y=1−e−x→e−x=1−y两边求对数可得:x=−ln(1−y)则 F − 1 ( x ) = − l n ( 1 − x ) F^{-1}(x)=-ln(1-x) F−1(x)=−ln(1−x),令 x i x_i xi为均匀分布样本,则 X i = − l n ( 1 − x i ) X_i=-ln(1-x_i) Xi=−ln(1−xi)为服从累积密度函数为F(x)分布的样本.
三、拒绝接受采样
拒绝接受采样的目的仍然是得到服从某个概率分布的样本,不过这种方法是直接利用概率密度函数(PDF)得到样本。如下图所示, p ( x ) p(x) p(x)是我们希望采样的分布, q ( x ) q(x) q(x)是我们提议的分布(proposal distribution), q ( x ) q(x) q(x)分布比较简单,令 k q ( x ) > p ( x ) kq(x)>p(x) kq(x)>p(x),我们首先在 k q ( x ) kq(x) kq(x)中按照直接采样的方法采样粒子,接下来判断这个粒子落在途中什么区域,对于落在蓝线以外的粒子予以拒绝,落在蓝线下的粒子接受,最终得到符合 p ( x ) p(x) p(x)的N个粒子。
拒绝接受采样的基本步骤:
- 生成服从q(x)的样本 → x i \rightarrow x_i →xi
- 生成服从均匀分布U(0,1)的样本—— u i u_i ui
- 当 q ( x i ) ⋅ u i < p ( x i ) q(x_i)\cdot u_i<p(x_i) q(xi)⋅ui<p(xi),也就是二维点落在蓝线以下,此时接受 X k = x i X_k=x_i Xk=xi
- 最终得到的 X k X_k Xk为服从p(x)的样本
实例
我们希望采样得到服从如下分布的样本
p ( x ) = e − x 2 2 ( sin 2 ( 6 + x ) + 3 cos 2 ( x ) sin 2 ( 4 x ) + 1 ) p(x) = e^{-\frac{x^2}{2}} \left( \sin^2(6+x) + 3 \cos^2(x) \sin^2(4x) + 1 \right) p(x)=e−2x2(sin2(6+x)+3cos2(x)sin2(4x)+1)
这个函数的形式比较复杂,不容易积分得到CDF的形式,也就更不好对CDF求反函数,所以我们采用拒绝-接受采样,直接对PDF进行采样.
定义拒绝抽样函数:
def Reject_Sampling(p,q_x):
"""对p(x)进行拒绝采样
Args:
p : 目标PDF函数
q_x : 选取的简单概率密度函数
Returns:
list: 采样样本
"""
size = 1e06 #初始采样点数目
k=10 #系数
sample=[]
for _ in range(int(size)): #不建议用for循环,速度慢,这样写比较好理解
xi=np.random.normal(0,1)
qi=k*q_x.pdf(xi) #乘上系数
ui=np.random.uniform(0,1)
pi=p(xi)
#判断
if qi*ui<pi:
sample.append(xi)
return sample
实验并作图:
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
def f(x):
px=np.exp(-(x**2)/2)*(np.sin(6+x)**2+3*(np.cos(x)**2)*(np.sin(4*x)**2)+1)
return px
q_x = stats.norm(0, 1)
RS=Reject_Sampling(f,q_x) #抽样数据集
fig,ax=plt.subplots(figsize=(10,5))
ax.hist(RS,bins=2000,density = True, stacked=True) #画出抽样分布
ax.set_xlabel('x', fontsize=15)
ax.set_ylabel("f(x)", fontsize=15)
ax.set_title("Rejection Sampling", fontsize=20)
我们发现采样的样本的分布图和目标函数分布图基本一致。
四、重要性采样
(1) 目的
重要性采样的目的:求一个f(x)在p(x)分布下的期望,即
E
[
f
(
x
)
]
=
∫
f
(
x
)
p
(
x
)
d
x
\mathbb{E}[f(x)]=\int f(x)p(x)dx
E[f(x)]=∫f(x)p(x)dx
当
p
(
x
)
p(x)
p(x)很复杂时,不解析,积分不好求时,可以通过重要性采样来计算。当
f
(
x
)
=
x
f(x)=x
f(x)=x,则可以算
E
[
x
]
\mathbb{E}[x]
E[x]在
p
(
x
)
p(x)
p(x)分布下的期望.
(2) 原理
首先, 当我们想要求一个函数 f ( x ) f(x) f(x) 在区间 [ a , b ] [a, b] [a,b] 上的积分 ∫ a b f ( x ) d x \int_{a}^{b} f(x) d x ∫abf(x)dx 时有可能会面临一个问题, 那就是积分曲线难以解析, 无法直接求积分。这时候我们可以采用一种估计的方式, 即在区 间 [ a , b ] [a, b] [a,b] 上进行采样: { x 1 , x 2 … , x n } \left\{x_{1}, x_{2} \ldots, x_{n}\right\} {x1,x2…,xn}, 值为 { f ( x 1 ) , f ( x 2 ) , … , f ( x n ) } \left\{f\left(x_{1}\right), f\left(x_{2}\right), \ldots, f\left(x_{n}\right)\right\} {f(x1),f(x2),…,f(xn)} .
如果采样是均匀的, 即如下图所示:
那么显然可以得到这样的估计: ∫ a b f ( x ) d x = b − a N ∑ i = 1 N f ( x i ) \int_{a}^{b} f(x) d x=\frac{b-a}{N} \sum_{i=1}^{N} f\left(x_{i}\right) ∫abf(x)dx=Nb−ai=1∑Nf(xi)在这里 b − a N \frac{b-a}{N} Nb−a 可以看作是上面小长方形的底部的 “宽”, 而 f ( x i ) f\left(x_{i}\right) f(xi) 则是坚直的 “长”。
上述的估计方法随着取样数的增长而越发精确,那么有什么方法能够在一定的抽样数量基础上来增加准确度,减少方差呢?比如x样本数量取10000,那么显然在f(x)比较大的地方,有更多的 x i x_i xi,近似的积分更精确,如下图所示,在圆圈部分 x i x_i xi样本应该更多,所以 x i x_i xi就不用均匀分布产生。
并且原函数 f ( x ) f(x) f(x) 也许本身就是定义在一个分布之上的, 我们定义这个分布为 p ( x ) p(x) p(x), 我们无法直接从 p ( x ) p(x) p(x) 上进行采样, 所以另辟蹊径重新找到一个更加简明的分布 q ( x ) q(x) q(x), 从它进行取样, 希望间接地求出 f ( x ) f(x) f(x) 在分布 p ( x ) p(x) p(x) 下的期望。
(2.1) 若 p ( x ) p(x) p(x)归一化
搞清楚了这一点我们可以继续分析了。首先我们知道函数
f
(
x
)
f(x)
f(x) 在概率分布
p
(
x
)
p(x)
p(x) 下的期望为:
E
[
f
(
x
)
]
=
∫
x
p
(
x
)
f
(
x
)
d
x
(
1
)
\mathbb{E}[f(x)]=\int_{x} p(x) f(x) d x \qquad(1)
E[f(x)]=∫xp(x)f(x)dx(1)
但是这个期望的值我们无法直接得到, 因此我们需要借助
q
(
x
)
q(x)
q(x) 来进行采样,
q
(
x
)
q(x)
q(x) 可以选取简单的分布,比如设q(x)为均匀分布,当我们在
q
(
x
)
q(x)
q(x) 上采样得到
{
x
1
,
x
2
,
…
,
x
n
}
\left\{x_{1}, x_{2}, \ldots, x_{n}\right\}
{x1,x2,…,xn} (即
x
i
x_i
xi服从q(x)分布)后,那么我们可以估计
f
f
f 在
q
(
x
)
\mathbf{q(x)}
q(x) 下的期望为:
E
[
f
(
x
)
]
=
∫
x
q
(
x
)
f
(
x
)
d
x
≈
1
N
∑
i
=
1
N
f
(
x
i
)
.
\mathbb{E}[f(x)]=\int_{x} q(x) f(x) d x \approx \frac{1}{N} \sum_{i=1}^{N} f\left(x_{i}\right) .
E[f(x)]=∫xq(x)f(x)dx≈N1i=1∑Nf(xi).
上面这个式子就简单很多了,只要我们得到
x
i
x_i
xi然后代入
f
(
x
)
f(x)
f(x)然后求和就行了,而且均匀分布的样本
x
i
x_i
xi很容易获得。接着我们来考虑原问题,对式(1)进行改写, 即:
p
(
x
)
f
(
x
)
=
q
(
x
)
p
(
x
)
q
(
x
)
f
(
x
)
p(x) f(x)=q(x) \frac{p(x)}{q(x)} f(x)
p(x)f(x)=q(x)q(x)p(x)f(x), 所以我们可以得到:
E
[
f
(
x
)
]
=
∫
x
q
(
x
)
p
(
x
)
q
(
x
)
f
(
x
)
d
x
\mathbb{E}[f(x)]=\int_{x} q(x) \frac{p(x)}{q(x)} f(x) d x
E[f(x)]=∫xq(x)q(x)p(x)f(x)dx
这个式子我们可以看作是函数
p
(
x
)
q
(
x
)
f
(
x
)
\frac{p(x)}{q(x)} f(x)
q(x)p(x)f(x) 定义在分布
q
(
x
)
q(x)
q(x) 上的期望, 当我们在
q
(
x
)
q(x)
q(x) 上采样
{
x
1
,
x
2
,
…
,
x
n
}
\left\{x_{1}, x_{2}, \ldots, x_{n}\right\}
{x1,x2,…,xn} (服从q(x)分布),可以估计
f
f
f 的期望:
E
[
f
(
x
)
]
=
1
N
∑
i
=
1
N
p
(
x
i
)
q
(
x
i
)
f
(
x
i
)
=
1
N
∑
i
=
1
N
w
i
f
(
x
i
)
\begin{aligned} \mathbb{E}[f(x)]&=\frac{1}{N} \sum_{i=1}^{N} \frac{p\left(x_{i}\right)}{q\left(x_{i}\right)} f\left(x_{i}\right)\\ &=\frac{1}{N} \sum_{i=1}^{N} w_i f\left(x_{i}\right) \end{aligned}
E[f(x)]=N1i=1∑Nq(xi)p(xi)f(xi)=N1i=1∑Nwif(xi)
在这里
w
i
=
p
(
x
i
)
q
(
x
i
)
w_i=\frac{p\left(x_{i}\right)}{q\left(x_{i}\right)}
wi=q(xi)p(xi)就是重要性权重。
(2.2)若 p ( x ) 没有归一化 p(x)没有归一化 p(x)没有归一化
上面的讨论是假设
p
(
x
)
p(x)
p(x)已经完成归一化了,也就是
∫
p
(
x
)
=
1
\int p(x)=1
∫p(x)=1,假如
p
(
x
)
p(x)
p(x)没有归一化,那么我们可以在上面的推导中对
p
(
x
)
p(x)
p(x)进行归一化:
E
[
f
(
x
)
]
=
∫
f
(
x
)
p
(
x
)
∫
p
(
x
)
d
x
d
x
=
∫
f
(
x
)
p
(
x
)
d
x
∫
p
(
x
)
d
x
=
∫
f
(
x
)
p
(
x
)
q
(
x
)
q
(
x
)
d
x
∫
p
(
x
)
q
(
x
)
q
(
x
)
d
x
.
\begin{aligned} \mathbb{E}[f(x)]&=\int f(x) \frac{p(x)}{\int p(x) d x} d x\\ &=\frac{\int f(x) p(x) d x}{\int p(x) d x}\\ &=\frac{\int f(x) \frac{p(x)}{q(x)} q(x) d x}{\int \frac{p(x)}{q(x)} q(x) d x}. \end{aligned}
E[f(x)]=∫f(x)∫p(x)dxp(x)dx=∫p(x)dx∫f(x)p(x)dx=∫q(x)p(x)q(x)dx∫f(x)q(x)p(x)q(x)dx.
而分子分母可分别得到,下面两式约等于都利用
q
(
x
)
q(x)
q(x)是均匀分布的假设:
∫
f
(
x
)
p
(
x
)
q
(
x
)
q
(
x
)
d
x
≈
1
n
∑
i
=
1
n
W
i
f
(
x
i
)
,
∫
p
(
x
)
q
(
x
)
q
(
x
)
d
x
≈
1
n
∑
i
=
1
n
W
i
.
\begin{aligned} \int f(x) \frac{p(x)}{q(x)} q(x) d x &\approx \frac{1}{n} \sum_{i=1}^{n} W_{i} f\left(x_{i}\right), \\ \int \frac{p(x)}{q(x)} q(x) d x &\approx \frac{1}{n} \sum_{i=1}^{n} W_{i}. \end{aligned}
∫f(x)q(x)p(x)q(x)dx∫q(x)p(x)q(x)dx≈n1i=1∑nWif(xi),≈n1i=1∑nWi.
其中
W
i
=
p
(
x
i
)
q
(
x
i
)
W_i=\frac{p(x_i)}{q(x_i)}
Wi=q(xi)p(xi),则最终可得
E
[
f
(
x
)
]
\mathbb{E}[f(x)]
E[f(x)]:
E
[
f
(
x
)
]
≈
∑
i
=
1
n
w
i
f
(
x
i
)
,
w
i
=
W
i
∑
i
=
1
n
W
i
\begin{aligned} \mathbb{E}[f(x)] \approx \sum_{i=1}^{n} w_{i} f\left(x_{i}\right), w_{i}=\frac{W_{i}}{\sum_{i=1}^{n} W_{i}} \end{aligned}
E[f(x)]≈i=1∑nwif(xi),wi=∑i=1nWiWi
(3)实例
我们希望求如下函数
p ( x ) = e − x 2 2 ( sin 2 ( 6 + x ) + 3 cos 2 ( x ) sin 2 ( 4 x ) + 1 ) , p(x) = e^{-\frac{x^2}{2}} \left( \sin^2(6+x) + 3 \cos^2(x) \sin^2(4x) + 1 \right) , p(x)=e−2x2(sin2(6+x)+3cos2(x)sin2(4x)+1),
的期望 E [ p ( x ) ] = ∫ − ∞ ∞ x p ( x ) d x \mathbb{E}[p(x)]=\int_{-\infty}^{\infty}xp(x)dx E[p(x)]=∫−∞∞xp(x)dx,这个函数的形式比较复杂,积分有点痛苦,所以我们尝试用重要性采样来得到这个期望值,注意在这个问题中 f ( x ) = x f(x)=x f(x)=x.
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
def p(x):
""" 概率密度函数,作为被积函数 """
px = np.exp(-(x**2)/2) * (np.sin(6 + x)**2 + 3 * (np.cos(x)**2) * (np.sin(4 * x)**2) + 1)
# px = 1 / 2 * np.exp(- (x-1) ** 2 /2) # 正太分布可以用来验证,均值是否为1
return px
def f(x):
return x
q_x = stats.norm(0, 1) # 这里我们取 q(x) 为标准正态分布
def importance_sampling(f,p, n=1000):
""" 使用重要性采样来估计函数的期望。
Args:
f (function): 目标函数 f(x)
p (function): f服从 p(x) 的分布,p(x)为概率密度函数
n (int): 采样的数量,默认为 1000。
Returns:
float: 估计的期望值 E[f(X)]。
"""
total_weighted_value = 0
total_weight = 0
for _ in range(n):
x_i = np.random.normal(0, 1) # 从 q(x) 中采样
W_i = p(x_i) / q_x.pdf(x_i) # 计算权重 W_i = p(x_i) / q(x_i)
total_weighted_value += f(x_i) * W_i
total_weight += W_i
expected_value = total_weighted_value / total_weight
return expected_value
estimated_mean = importance_sampling(f, p, 100000)
print("E(X) =", estimated_mean)
最终我们得到期望值为
E
[
p
(
x
)
]
=
−
0.0343
\mathbb{E}[p(x)]=-0.0343
E[p(x)]=−0.0343,我们还可以用上面拒绝接受采样的结果验证,拒绝接受采样是得到服从
p
(
x
)
p(x)
p(x)分布的样本,那么我们可以用样本均值来估计期望,即:
μ
=
x
1
+
x
2
+
⋯
+
x
n
n
\mu=\frac{x_1+x_2+\cdots+x_n}{n}
μ=nx1+x2+⋯+xn
mean = np.mean(RS)
print("mean:",mean)
我们得到 μ = − 0.03178 \mu=-0.03178 μ=−0.03178,可以发现两个的结果一致(采样有随机性,每次结果都有细微差别)。我们还可以用正太分布来验证,比如代码里我写的 p ( x ) = 1 2 e − ( x − 1 ) 2 2 p(x)=\frac12e^{-\frac{(x-1)^2}{2}} p(x)=21e−2(x−1)2,显然这个正太分布的期望为1,我们用重要性采样得到的结果为 μ = 1.0020 \mu=1.0020 μ=1.0020,可以看到基本上可以精确到小数点后三位,相对误差为0.206%.