取樣演算法簡述

1 閱讀之前

在閱讀本文之前請先明確什麼是樣本。

前置：假設一個樣本取自分佈p（x），那麼我們通常將其記錄為

$X\sim p(x)$

，但本文使用的小寫字母加下標的方式表示樣本

$x_i\sim p(x)$

。

前置：單個樣本是無法研究分佈的，只有樣本數量足夠時才能研究性質。

2 均勻分佈隨機數的產生

在獲取均勻分佈的隨機數後可以通過後續的演算法來獲取其他分部的隨機數。最簡單的方式是對座標進行變換。假設

是原始空間中的座標，而

是變換空間後的座標。定義函式

$f(x),x\in \mathbb{R}^2$

和座標變換：

此時假設函式

為機率密度函式，而到目前為止僅能獲取均勻分佈的隨機數，因此機率密度函式定義為

定義變換：

$\left\{ \begin{aligned} y_1&=&\sqrt{-2\ln{(x_1)}}\cos(2\pi x_2)\\ y_2&=&\sqrt{-2\ln{(x_2)}}\sin(2\pi x_1) \end{aligned} \right.$

機率密度函式有積分意義：

$F=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}fds$

變換域面積與原空間面積關係為：

其中

$J_{ij}=\frac{\partial y_i}{\partial x_i}$

為雅克比，因此函式可以轉換為：

$F=\int f(x)dx_1 dx_2=\int f(y) |J|dy_1 dy_2=\int \frac{1}{\sqrt{2\pi}}\exp{(-y_1^2/2)}\frac{1}{\sqrt{2\pi}}\exp{(-y_1^2/2)} dy_1 dy_2$

所以均勻分佈的樣本可以產生標準正態分佈的樣本。

對於座標變換而言，還有更加通用的方式。假設

是（0-1）區間均勻分佈的樣本。而我們想獲取的分佈為

。記錄：

$F(y)=\int_{-\infty}^yp(\theta)d\theta$

定義座標變換

變數

的累計分佈函式為：

$F(x)=\int_0^x 1 d\theta=\int_{-\infty}^y \frac{\partial x}{\partial y} d\phi=\int_{-\infty}^y p(\phi) d\phi$

所以如果

$y=F^{-1}(x)$

的話就是所需分佈樣本。對於難以計算的情況，我們可以使用近似的方式來計算機率：

import

numpy

def

（

）：

“”“

所需機率密度函式

”“”

return

。

exp

（

）

。

def

sampler

（）：

“”“

生成一個隨機樣本

”“”

# 產生多個樣本

。

random

。

random

（［

］）

# 計算機率

（

）

# 計算累計機率

。

cumsum

（

）

。

sum

（

）

# 產生均勻分佈隨機數

。

random

。

random

（）

# 根據機率選擇樣本

sample_id

int

（

。

searchsorted

（

，

。

random

。

random

（）））

return

［

sample_id

］

# 給定初始值

x_samples

［］

n_samples

5000

for

itr

range

（

n_samples

）：

x_samples

。

append

（

sampler

（））

x_samples

。

array

（

x_samples

）

import

matplotlib。pyplot

plt

。

hist

（

x_samples

）

plt

。

show

（）

所得結果

3 透過拒絕部分樣本獲取所需分佈

3.1 拒絕接受取樣

均勻分佈樣本與其他分佈樣本區別就在於一些低機率樣本出現過多，因此可以使用拒絕部分樣本的方式對樣本進行取樣。假設需要的分佈為

。其他分部樣本為

。二者之間關係為

$p(x)=cq(x)\alpha(x)$

為使得

在之前乘以一個乘子。那麼取樣過程為：

從分佈q中產生樣本

根據接受率

$\alpha(x)=\frac{p(x)}{cq(x)}$

接受樣本

產生均勻分佈隨機數

$s\sim U(0,1)$

，如果

$s<\alpha(x)$

則接受樣本。

根據這個方法舉個例子，產生二維正態分佈的隨機數。

“”“

用於演示拒絕接受取樣演算法

”“”

import numpy as np

def p（x）：

“”“

所需機率密度函式

”“”

return np。exp（-np。sum（x**2）/2）/np。pi/2

def q（x）：

“”“

原始樣本的機率密度函式

”“”

return 1/100

# 產生均勻分佈的樣本

x = np。random。random（［10000， 2］） * 10 - 5

# 選擇後樣本

x_accept = ［］

x_reject = ［］

c = 100/np。pi/2

for itr in x：

alpha = p（itr） / q（itr） / c

# 產生閾值

s = np。random。random（）

if s < alpha：

x_accept。append（itr）

else：

x_reject。append（itr）

x_accept = np。array（x_accept）

x_reject = np。array（x_reject）

# 繪圖

import matplotlib。pyplot as plt

plt。subplot（121）

x_norm = np。random。normal（0， 1，［1000， 2］）

plt。scatter（x_norm［：， 0］， x_norm［：， 1］）

plt。xlim（［-5， 5］）

plt。ylim（［-5， 5］）

plt。subplot（122）

plt。scatter（x_accept［：， 0］， x_accept［：， 1］， marker=‘o’， color=“#000099”， alpha=1。0）

plt。scatter（x_reject［：， 0］， x_reject［：， 1］， marker=‘+’， color=“#990000”， alpha=0。5）

plt。xlim（［-5， 5］）

plt。ylim（［-5， 5］）

plt。show（）

結果

3.2 Metropolis-Hastings取樣

透過執行上述程式可以看到，透過拒絕接受樣本在計算過程中有大量計算是無用的。因此選取合適的分佈

是十分必要的。但這通常都很難獲得。由此產生出了新的Metropolis-Hastings（MH）演算法。下面對演算法進行詳細描述。

在拒絕接受取樣的過程中發現，每個樣本產生的過程是完全獨立的。這符合隨機數產生的邏輯。但如果將取樣看作一個過程，那麼隨機數產生的過程並非完全獨立的。。依然以上面的例子來說，我們產生了一個均勻分佈的隨機數。而我們所期望的分佈是標準的二維正態分佈。這需要我們拒絕部分樣本。到此為止與拒絕接受取樣都沒有什麼不同。但是假設我們已經產生了一個機率非常小的樣本x，那麼之後的樣本應當比x大一些才能平衡採到x所帶來的影響。由此可以將取樣看作一個過程：

$\begin{aligned} x_t & \sim & p_t(x)\\ x_{t+1} & \sim & p_{t+1}(x) \end{aligned}$

作為取樣過程，我們希望在取樣過程中機率不變，也就是一個平穩的過程：

$p_{t+1}(x)=p_t(x)=p(x)$

這裡定義一個轉移矩陣

$Q(x_t,x_{t+1})=p(x_{t+1}|x_t)$

：

$p(x_{t+1})=\int_{-\infty}^{+\infty}p(x_t)Q(x_t,x_{t+1})dx_t$

當

$p(x_{t+1})Q(x_{t+1},x_t)=p(x_t)Q(x_t, x_{t+1})$

時，上式變為：

$p(x_{t+1})=\int_{-\infty}^{+\infty}p(x_{t+1})p(x_t| x_{t+1})dx_t=p(x_{t+1})$

等式恆成立，因此條件也稱為細緻平穩條件。滿足細緻平穩條件的樣本，機率是分佈是平穩的。這個過程並不是完整的證明。通常轉移機率

是非常難以獲得的。因此我們會構造一個符合細緻平穩條件的

，比如

$p(x_{t+1})Q(x_{t+1},x_t)\frac{p(x_t)Q(x_t, x_{t+1})}{\max{(p(x_t)Q(x_t, x_{t+1}),p(x_{t+1})Q(x_{t+1},x_t))}}=p(x_t)Q(x_t, x_{t+1})\frac{p(x_{t+1})Q(x_{t+1},x_t)}{\max{(p(x_t)Q(x_t, x_{t+1}),p(x_{t+1})Q(x_{t+1},x_t))}}$

定義

$\alpha=\frac{p(x_{t+1})Q(x_{t+1},x_t)}{\max{(p(x_t)Q(x_t, x_{t+1}),p(x_{t+1})Q(x_{t+1},x_t))}}=min(\frac{p(x_{t+1})Q(x_{t+1},x_t)}{p(x_{t})Q(x_{t},x_{t+1})}, 1)$

因此Metropolis-Hastings演算法為

從條件機率

$Q(x,x_{t+1})$

中產生新樣本

$x_{t+1}$

產生均勻分佈

$s\sim U(0,1)$

如果

$s<\alpha$

則接受取樣

使用上述演算法，產生二維正態分佈的隨機數。原始樣本產生過程中我們依然使用均勻分佈樣本。這些樣本每個採到的機率均是相同的，因此轉移矩陣

$Q(x_{t+1},x_{t})=Q(x_{t},x_{t+1})$

此時

$\alpha=\min(\frac{p(x_{t+1})}{p(x_t)},1)$

“”“

用於演示Metropolis-Hastings演算法

”“”

import numpy as np

def p（x）：

“”“

所需機率密度函式

”“”

return np。exp（-np。sum（x**2）/2）/np。pi/2

def cal_alpha（x1， x2）：

return np。min（［p（x2）/p（x1）， 1］）

# 產生均勻分佈的樣本

x = np。random。random（［1000， 2］） * 10 - 5

# 選擇後樣本

x_accept = ［］

x_reject = ［］

c = 100/np。pi/2

x_accept。append（x［0］）

for itr in range（len（x）-1）：

alpha = cal_alpha（x_accept［-1］， x［itr+1］）

s = np。random。random（）

if s < alpha：

x_accept。append（x［itr+1］）

else：

x_accept。append（x_accept［-1］）

x_reject。append（x［itr+1］）

x_accept = np。array（x_accept）

x_reject = np。array（x_reject）

# 繪圖

import matplotlib。pyplot as plt

plt。subplot（121）

x_norm = np。random。normal（0， 1，［1000， 2］）

plt。scatter（x_norm［：， 0］， x_norm［：， 1］）

plt。xlim（［-5， 5］）

plt。ylim（［-5， 5］）

plt。subplot（122）

plt。scatter（x_accept［：， 0］， x_accept［：， 1］， marker=‘o’， color=“#000099”， alpha=1。0）

plt。scatter（x_reject［：， 0］， x_reject［：， 1］， marker=‘+’， color=“#990000”， alpha=0。5）

plt。plot（x_accept［：， 0］， x_accept［：， 1］， color=“#000099”， alpha=0。3）

plt。xlim（［-5， 5］）

plt。ylim（［-5， 5］）

plt。show（）

計算結果

3.3 吉布斯取樣

上面的取樣過程與狀態轉移的方式以及轉移矩陣的選擇上實際上並無太大關係。可以發現MH演算法的的基本過程就是選擇一個樣本

並計算其機率

，如果新生成一個樣本

其機率

，如果新的樣本j的機率大於i的機率那麼我們就接受它，因為比它機率小的樣本都接受了。如果j的機率小於i的機率，那麼我們就按照二者機率之差來隨機決定是否接受此樣本。這個過程依然會拒絕樣本，在高維時則更加難以取樣。

現在我們假設取樣過程中先後生成了兩個樣本

$x_t,x_{t+1}$

，二者僅在座標i上是不同的。此時二者機率為

$\left\{ \begin{aligned} p(x_t)&=&p(x_t^{\neg i})p(x_t^i|x_t^{\neg i})\\ p(x_{t+1})&=&p(x_{t+1}^{\neg i})p(x_{t+1}^i|x_{t+1}^{\neg i}) \end{aligned} \right.$

其中

$x_t^{\neg i}$

代表x中不含i的元素，由於只有i不同，所以

$x_t^{\neg i}=x_{t+1}^{\neg i}=x^{\neg i}=a$

。上式是我們對聯合機率進行的分解。由於

$p(x_{t+1}^{\neg i})=p(x_{t}^{\neg i})$

，所以

$p(x_t)p(x^i_{t+1}|x^{\neg i})=p(x_{t+1})p(x^i_t|x^{\neg i})$

這是符合細緻平穩條件的。

“”“

用於演示吉布斯取樣演算法

”“”

import numpy as np

def p（x）：

“”“

所需機率密度函式

”“”

return np。exp（-np。sum（x**2， axis=1）/2）/np。pi/2

def dim_sampler（x， dim）：

“”“

從某個維度取樣

直接計算多個樣本的機率，根據機率選擇合適的樣本。

由於僅是一維分佈，我們也可以用其他方式進行取樣。

”“”

# 產生多個樣本

sp = np。random。random（［20］） * 10 - 5

xs = np。ones（［20， len（x）］）*np。array（x）

xs［：， dim］ = sp

# 計算機率

ps = p（xs）

# 根據機率選擇樣本

t = np。cumsum（ps）/np。sum（ps）

s = np。random。random（）

sample_id = np。searchsorted（t， s）

return list（xs［sample_id］）

# 給定初始值

x_samples = ［［0， 0］］

n_samples = 600

state = x_samples［-1］

for itr in range（n_samples）：

for dim in range（len（state））：

state = dim_sampler（state， dim）

x_samples。append（state）

x_samples = np。array（x_samples）

import matplotlib。pyplot as plt

plt。subplot（121）

x_norm = np。random。normal（0， 1，［1000， 2］）

plt。scatter（x_norm［：， 0］， x_norm［：， 1］）

plt。xlim（［-5， 5］）

plt。ylim（［-5， 5］）

plt。subplot（122）

plt。scatter（x_samples［：， 0］， x_samples［：， 1］， marker=‘o’， color=“#000099”， alpha=1。0）

plt。plot（x_samples［：， 0］， x_samples［：， 1］， color=“#000099”， alpha=0。3）

plt。xlim（［-5， 5］）

plt。ylim（［-5， 5］）

plt。show（）

所得結果

對於吉布斯取樣來說，如果有多個峰，那麼取樣過程中可能出現轉移機率過低，無法採集到的問題：

小蜜蜂問答

小蜜蜂問答

取樣演算法簡述

推薦文章

小蜜蜂問答

小蜜蜂問答

取樣演算法簡述

相關文章

農作物產量的調查方法包括哪兩大類？

knn演算法的原理與實現

利用python進行相關性分析

如何用最少的操作打出 N 個 1？

推薦文章