随機采樣
采樣是根據某種分布去生成一些資料點。最基本的假設是認為我們可以獲得服從均勻分布的随機數，再根據均勻分布生成複雜分布的采樣。對于離散分布的采樣，可以把機率分布向量看作一個區間段，然後判斷u落在哪個區間段内。對于比較複雜的分布比如正态分布我們可以通過Box-Muller算法，實作對高斯分布的采樣。
其他幾個著名的連續分布，包括指數分布，Gamma分布，t分布，F分布，Beta分布，Dirichlet分布都可以通過類似的數學變換得到。
蒙特卡洛估計
從一個分布中采樣可以做很多有用的事情，包括求邊緣分布，MAP推斷，以及計算如下形式的期望
拒絕采樣
舉例：比如要計算圓的面積可以在矩形中均勻的撒點，統計落在圓中的點的個數處于總的點的個數，就是圓的面積與矩形面積的比值。
假設我們已經可以抽樣高斯分布q(x)（如Box–Muller_transform 算法），我們按照一定的方法拒絕某些樣本，達到接近p(x)分布的目的。
拒絕采樣時蒙特卡洛中的一種特殊情況。計算p(x=x'):我們可以将這個機率寫成如下形式
接着采用蒙特卡洛近似
重要采樣
拒絕采樣會浪費大量的樣本，一個更好的方法是通過重要采樣進行采樣。
本來需要從p分布中進行采樣，由于種種原因（p分布不好采樣等）我們不願意從p分布中采樣，而是從q分布中采樣，需要在f中乘上權重因子w(x)
對于p(x=x')可以通過如下的計算
Python代碼
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
import seaborn
def qsample():
return np.random.rand() * 4.
def p(x):
return 0.3 * np.exp(-(x - 0.3) ** 2) + 0.7 * np.exp(-(x - 2.) ** 2 / 0.3)
def q(x):
return 4.0
def importance(nsamples):
samples = np.zeros(nsamples, dtype=float)
w = np.zeros(nsamples, dtype=float)
for i in range(nsamples):
samples[i] = qsample()
w[i] = p(samples[i]) / q(samples[i])
return samples, w
def sample_discrete(vec):
u = np.random.rand()
start = 0
for i, num in enumerate(vec):
if u > start:
start += num
else:
return i - 1
return i
def importance_sampling(nsamples):
samples, w = importance(nsamples)
# print samples
final_samples = np.zeros(nsamples, dtype=float)
w = w / w.sum()
# print w
for j in range(nsamples):
final_samples[j] = samples[sample_discrete(w)]
return final_samples
x = np.arange(0, 4, 0.01)
x2 = np.arange(-0.5, 4.5, 0.1)
realdata = 0.3 * np.exp(-(x - 0.3) ** 2) + 0.7 * np.exp(-(x - 2.) ** 2 / 0.3)
box = np.ones(len(x2)) * 0.8
box[:5] = 0
box[-5:] = 0
plt.plot(x, realdata, 'g', lw=6)
plt.plot(x2, box, 'r--', lw=6)
# samples,w = importance(5000)
final_samples = importance_sampling(5000)
plt.hist(final_samples, normed=1, fc='c')
plt.show()
MCMC（Markov Chain Monte Carlo）
MCMC采樣算法
從馬爾科夫鍊中對某個分布進行采樣，采樣的結果通過蒙特卡洛方法來對期望進行估計。
一個直覺的了解是：收集了可能來自需要觀測山的大量地點，我們就可以近似的構造出這座山出來。
從目前地點開始，移動到一個靠近你的新地點。如果該地點可能來自需要觀測的山，那麼移動到該地點，如果不是待在原處。通過大量的疊代，傳回所有可能的地點。
MCMC适合使用在高維的空間中，會提高采樣的效率。
細緻平穩條件
假設我們已經有一個轉移矩陣為Q馬氏鍊（q(i,j)）表示從狀态i轉移到狀态j的機率，通常情況下
也就是細緻平穩條件不成立，是以p(x)不太可能是這個馬氏鍊的平穩分布。為此需要引入
是以有
Metropolis-Hastings采樣算法
以上MCMC算法已經可以工作了，但是它有一個小的問題：馬氏鍊Q在轉移的過程中可能偏小，這樣采樣過程中馬氏鍊容易原地踏步。是以需要提高接受率。
假設
此時滿足細緻平穩條件，于是
上式兩邊擴大5倍，我們改寫為
我們提高了接收率，而細緻平穩條件沒有打破。這啟發我們可以把細緻平穩條件的同比例放大，使得兩數中最大的一個放大到1，這樣就提高了采樣中的跳轉接受率。
Python代碼：
# The Metropolis-Hastings algorithm
from pylab import *
from numpy import *
def p(x):
mu1 = 3
mu2 = 10
v1 = 10
v2 = 3
return 0.3*exp(-(x-mu1)**2/v1) + 0.7* exp(-(x-mu2)**2/v2)
# return -x**2
def q(x):
mu = 5
sigma = 10
return exp(-(x-mu)**2/(sigma**2))
stepsize = 0.5
x = arange(-10,20,stepsize)
px = zeros(shape(x))
for i in range(len(x)):
px[i] = p(x[i])
N = 5000
# independence chain
u = random.rand(N)
mu = 5
sigma = 10
y = zeros(N)
y[0] = random.normal(mu,sigma)
for i in range(N-1):
ynew = random.normal(mu,sigma)
alpha = min(1,p(ynew)*q(y[i])/(p(y[i])*q(ynew)))
if u[i] < alpha:
y[i+1] = ynew
else:
y[i+1] = y[i]
# random walk chain
u2 = random.rand(N)
sigma = 10
y2 = zeros(N)
y2[0] = random.normal(0,sigma)
for i in range(N-1):
y2new = y2[i] + random.normal(0,sigma)
alpha = min(1,p(y2new)/p(y2[i]))
if u2[i] < alpha:
y2[i+1] = y2new
else:
y2[i+1] = y2[i]
figure(1)
nbins = 30
hist(y, bins = x)
plot(x, px*N/sum(px), color='g', linewidth=2)
plot(x, q(x)*N/sum(px), color='r', linewidth=2)
figure(2)
nbins = 30
hist(y2, bins = x)
plot(x, px*N/sum(px), color='g', linewidth=2)
show()
該代碼用到一個比較特殊的馬爾科夫過程：當一個随機變量，其增量的變化服從于正态分布時，我們說，這個随機過程叫做維納過程（Wiener process)——或者一個更出名的名字：布朗運動過程
設x(t)是一個獨立增量過程，即随機過程x(t)的下一階段增量的機率分布獨立于其他任意階段的機率分布
X(t)是一個關于t的連續函數；
X(t)的增量服從關于時間t的正态分布，即
Gibbs采樣
由于Metropolis-Hastings采樣的通常小于1。是以Gibbs進一步提高到1。
考慮二維的情形，假設有一個機率分布p(x,y)，考慮x坐标相同的兩個點A(x1,y1),B(x1,y2)
我們發現
是以得到
即
基于以上等式，我們發現，在x=x1這條平行于y軸的直線上，如果使用條件分布p(y|x1)做為任何兩點之間的轉移機率，那麼任何兩個點之間的轉移滿足細緻平穩條件。
于是我們可以如下構造平面上任意兩點之間的狀态轉移機率矩陣Q
有了如上的轉移矩陣Q，我們很容易驗證對平面上任意兩點X，Y，滿足細緻平穩條件
以上采樣過程中，如圖所示，馬氏鍊的轉移隻是輪換的沿着坐标軸x軸和y軸做轉移
很容易将算法擴充到高維的情況中
def pXgivenY(y,m1,m2,s1,s2):
return np.random.normal(m1 + (y-m2)/s2,s1)
def pYgivenX(x,m1,m2,s1,s2):
return np.random.normal(m2 + (x-m1)/s1,s2)
def gibbs(N=5000):
k=20
x0 = np.zeros(N,dtype=float)
m1 = 10
m2 = 20
s1 = 2
s2 = 3
for i in range(N):
y = np.random.rand(1)
# 每次采樣需要疊代 k 次
for j in range(k):
x = pXgivenY(y,m1,m2,s1,s2)
y = pYgivenX(x,m1,m2,s1,s2)
x0[i] = x
return x0
def f(x):
"""目标分布"""
return np.exp(-(x-10)**2/10)
# 畫圖
N=10000
s=gibbs(N)
x1 = np.arange(0,17,1)
plt.hist(s,bins=x1,fc='c')
x1 = np.ar
ange(0,17, 0.1)
px1 = np.zeros(len(x1))
for i in range(len(x1)):
px1[i] = f(x1[i])
plt.plot(x1, px1*N*10/sum(px1), color='r',linewidth=3)
plt.show()