馬爾可夫性
過程或者系統在時刻t0處的狀態為已知條件下，過程(或系統)在時刻t>t0所處狀態的條件分布與過程在時刻t0之前所處的狀態無關(即，在已知過程"現在"的條件下，其"將來"不依賴於"過去")
馬爾可夫性的函數表示
1.隨機過程：設T是以無限實數集，我們把依賴於參數t∈T的一族(無限多個)隨機變量成為隨機過程，記為{X(t), t∈T}
對於每個t∈T，X(t)是一個隨機變量，T叫做參數集合
若將t看作為時間，稱X(t)為：時刻t時過程的狀態，X(ti)=x(實數)稱為t=ti時，過程處於狀態x，對於一切t∈T，X(t)所有可能取得一切值的全體成為隨機過程的狀態空間
2.函數表示
設隨機過程{X(t), t∈T}的撞塌其空間為I，如果對時間t的任意n個數值t1<t2<.<tn, n >= 3, ti ∈ T, 在條件X(ti)=xi, xi ∈ I, i=1~n-1下，X(tn)的條件分布函數恰等於在條件X(tn-1)=xn-1下X(tn)的條件分布函數
即：
P{X(tn)<=xn | X(t1)=x1, X(t2)=x2, . ,X(tn-1)=xn-1}
=> P{X(tn)<=xn | X(tn-1)=xn-1}, xn ∈ R
或者寫為：
Ftn|t1.tn-1(Xn, tn | x1, x2, . , xn-1; t1, t2, ., tn-1)
= Ftn|tn-1(xn, tn | xn-1, tn-1)
該過程{X(t), t∈T}具有馬爾可夫性和成為無後效性，並稱此過程為馬爾可夫過程
時間和狀態都是離散的馬爾可夫過程稱為馬爾可夫鏈（Markov process --> Markov chain），簡稱為馬氏鏈
馬爾可夫鏈記為{Xn=X(n), n=0, 1, 2, .}，它可以理解為看作在時間集合T1={0, 1, 2, .}上對離散狀態的馬爾可夫過程進行時序上相繼觀察的結果，約定記鏈的狀態空間為I={a1, a2, . }, ai ∈ R，在鏈的情形，馬爾可夫性通常用條件分布律來表示
即對任意的正整數n, r和0<=t1<t2<.<tr<m, n+m ∈ T1, 有：
P{Xm+n = aj | Xt1 = ai1, Xt2 = ai2, ., Xtr = air, Xm = ai}
= P{Xm+n = aj | Xm = ai}, 其中 a. ∈ I.
記上式右端為Pij(m, m+n)，成條件概率Pij(m, m+n) = P{Xm+n = aj | Xm = ai}為馬爾可夫鏈在時刻m處處於狀態ai的條件下，在時刻m+n轉移到狀態aj的轉移概率
由於馬爾可夫鏈在時刻m從任何一個狀態ai出發到另一個時刻m+n，必然轉移到a1, a2, . 狀態中的某一個，所以 ∑(j=1~+∞)Pij(m, m+n)=1, i=1, 2, .
由轉移概率組成的矩陣P(m, m+n)=(Pij(m, m+n))稱為馬爾可夫鏈的轉移概率矩陣，該矩陣每一行的元之和等於1
當轉移概率Pij(m, m+n)只與i,j及時間間距n有關時，可以將其記為Pij(n),即
Pij(m, m+n) = Pij(n) [*只看重轉移的時間間距n]
並稱該轉移概率具有平穩性，同時也稱該鏈是齊次的或時齊的
*馬爾可夫鏈為齊次的情形下，之前定義的轉移概率：
Pij(m, m+n) = Pij(n) = P{Xm+n = aj | Xm = ai}
則上述轉移概率稱為馬爾可夫鏈的n步轉移概率，P(n) = (Pij(n))為n步轉移概率矩陣
*一步轉移概率矩陣：Pij=Pij(1)=P{Xm+1=aj | Xm = ai}
*一步轉移概率矩陣：
Xm\Xm+1 a1 a2 aj ******
a1 [P11 P12 ****** P1j ******]
a2 [P21 P22 ****** P2j ******]
* [*** *** ****** *** ******] = P(1) ==>(記作)P.
* [*** *** ****** *** ******]
ai [Pi1 Pi2 ****** Pij ******]
狀態轉移矩陣的行(Xm+1的狀態到aj) 狀態轉移矩陣的列(Xm的狀態到ai)
上述矩陣左側和上側標記上狀態a1, a2, . 為了顯示Pij是由狀態ai經過一步轉移至aj的概率
馬爾可夫鏈的有限維分布
首先，記pj(0) = P{X0 = aj}, aj ∈ I. j=1,2,. (0時刻時，系統處於aj狀態下的概率) 稱其為馬爾可夫鏈的初始分布
再看馬爾可夫鏈在任一時刻n∈T1的一維分布 pj(n) = P{Xn = aj}, aj ∈ I, j=1,2,. (n時刻時，系統處於aj狀態下的概率)
顯然，應該有∑(j=1~+∞)pj(n) = 1 (n時刻時，系統處於a1~a∞狀態下的概率累加為1)
且又有：P{Xn = aj} = ∑(i=1~+∞)P(X0 = ai, Xn = aj)
= ∑(i=1~+∞)P{Xn = aj | X0 = ai}*P{X0 = ai}
或者，pj(n) = ∑(i=1~+∞)pi(0)Pij(n), j=1,2,3, .
馬爾可夫鏈的一維分布也可以用行向量形式表示：
P(n) = [P{Xn = a1}, P{Xn = a2}, ., P{Xn = aj}] => [p1(n), p2(n), ., pj(n), . ]
利用矩陣乘法(I為可列無限集合時，仍用有限階矩陣乘法的規則確定矩陣的元)
則有：pj(n) = ∑(i=1~+∞)pi(0)*Pij(n) => p(n) = p(0)P(n)
*P(n) = (Pij(n)) = Pij(n) = P{Xm+n = aj | Xm = ai} [P(n)為n步轉移概率矩陣]
*上式表明，馬爾可夫鏈在任一一時刻n∈T1時的一維分布由初始分布p(0)和n步轉移的概率矩陣所確定
*且又存在，對於任意n個時刻t1<t2< .. <tn, ti ∈ T1，以及狀態ai1, ai2, ., ain ∈ I，馬爾可夫鏈的n維分布：
P{Xt1 = ai1, Xt2 = ai2, ., Xtn = ain} = P{Xt1 = ai1}P{Xt2=ai2 | xt1 = ai1} * . *P{Xtn = ain | Xt1 = ai1, Xt2 = ai2, ., Xtn-1 = ain-1}
= P{Xt1 = ai1}P{Xt2 = ai2 | Xt1 = ai1}* . *P{Xtn = ain | Xtn-1 = ain-1}
=> pi1(t1)Pi1i2(t2 - t1) . Pin-1in(tn - tn-1)
*可知：馬爾可夫鏈的有限維分布同樣完全由初始分布和轉移概率所決定(轉移概率決定了馬爾可夫鏈運動的統計規律)
多步轉移概率確定
馬爾可夫鏈的n步轉移概率Pij(n)滿足的基本方程
設{X(n), n=0,1,2, .}是一個齊次馬爾可夫鏈，則對任意的u, v∈T1，有
Pij(u+v) = ∑(k=1~+∞)Pik(u)Pkj(v), i,j = 1,2,3, .
上式即為「查普曼-科爾莫戈羅夫」方程(或簡稱C-K方程)
*C-K方程：即「從時刻s所處的狀態ai，即X(s)=ai出發，經過時間段u+v轉移到狀態aj，即X(s+u+v)=aj」，該事件可分解為「從X(s)=ai出發，先經過時間段u轉移至中間狀態ak(k=1,2,3,.)，再從ak經時間段v轉移到狀態aj」 這樣一系列事件的和事件
*Chapman-Kolmogorov(查普曼-科爾莫戈羅夫)方程的證明過程
證明：先固定ak∈I, 和s∈T1，由條件概率定義和乘法定理，有
P{X(s+u+v) = aj, X(s+u) = ak | X(s) = ai}
=> P{X(s+u) = ak | X(s) = ai}P{X(s+u+v) = aj | X(s+u) = ak, X(s) = ai}
=> Pik(u)Pkj(v) [馬爾可夫性和齊次性遵循！]
hint:Pij(m, m+n) = Pij(n) = P{Xm+n = aj | Xm = ai}
且又由於事件組「X(s+u) = ak」, k=1,2,3, . 構成一個劃分，故存在
Pij(u+v) = P{X(s+u+v) = aj | X(s) = ai}
= ∑(k=1~+∞)P{X(s+u+v) = aj , X(s+u) = ak | X(s) = ai}
將P{X(s+u+v) = aj, X(s+u) = ak | X(s) = ai} = Pik(u)Pkj(v)代入上式中得
Pij(u+v) = ∑(k=1~+∞)Pik(u)Pkj(v) 方程得證
*對Chapman-Kolmogorov(查普曼-科爾莫戈羅夫)方程的圖例補充
*C-K方程的矩陣形式
P(u+v) = P(u)P(v) P(*)為*步概率轉移矩陣
可以利用C-K方程易確認n步轉移概率，在上式中令u=1, v=n-1，得到遞推關係
P(n) = P(1)P(n-1) = PP(n-1)
-->P(n-1) = P(1)P(n-2) = PP(n-2)
-->依此類推 得P(n) = P^n
*則可以確定，對齊次馬爾可夫鏈而言，n步轉移概率矩陣是一步轉移概率矩陣的n次方
*只有兩個狀態的馬爾可夫鏈，一步轉移概率矩陣一般可以表示為：
0 1
P = 0 [1-a a]
1 [b 1-b] 0< a,b <1 0 1
可得n步轉移概率矩陣為P(n) = P^n = P = 0 [P00(n) P01(n)]
1 [P10(n) P11(n)]
==>1/(a+b)[b a] + (1-a-b)^n/(a+b)[a -a]
[b a] [-b b] n = 1,2,3, .
馬爾可夫鏈的遍歷性
對於一般的2個狀態的馬爾可夫鏈，由之前上式得知，當0<a,b<1時，Pij(n)有極限
lim(n->+∞) P00(n) = lim(n->+∞) P10(n) = b/(a+b) ==>記為π0
lim(n->+∞) P01(n) = lim(n->+∞) P11(n) = a/(a+b) ==>記為π1
*上述極限的定義：對固定的狀態j，不管鏈在某一時刻從什麼狀態(i=0或1)出發，通過長時間的狀態轉移，到達狀態j的概率都趨近於πj，即所謂的遍歷性
又由於π0+π1=1，故(π0, π1) ==(記作)>π 構成一個分布律，稱其為鏈的極限分布
*一般地，設齊次馬爾可夫鏈的狀態空間為I，若對於所有的ai,aj ∈ I， 轉移概率Pij(n)存在極限
lim(n->+∞)Pij(n) = πj (不依賴於i)
或者用矩陣形式表示為：P(n) = P^n ---->(n->+∞)
[π1, π2, π3, ., πj, .]
[π1, π2, π3, ., πj, .]
[π1, π2, π3, ., πj, .]
[...]
則稱此鏈具有遍歷性，又若∑(j)πj = 1，則同時稱 π = [π1, π2, . ]為鏈的極限分布
*定理：設齊次馬爾可夫鏈{Xn, n>=1}的狀態空間為I={a1,a2,a3, ., aN}，P是他的一步轉移概率矩陣，如果存在正整數m，便對任意的ai, aj ∈ I，都有
Pij(m) > 0, i,j = 1,2,3, . , N
則此鏈具有遍歷性，且有極限分布π = [π1, π2, . , πN]它是方程組π = πP 或 πj = ∑(j=1~N)πiPij, j = 1,2,3, . ,N的滿足條件 πj > 0, ∑(j=1~N)πj=1的唯一解
上述定理證明，證明有限鏈是遍歷的，只需要找到一個正整數m，使m步轉移概率矩陣P^m無零元，在該定理的條件下，馬爾可夫鏈的極限分布又是平穩分布，即p(0) = π，則鏈在任一時刻n∈T1的分布p(n)永遠與π一致
*由之前定理：
p(n) = p(0)P(n) P(n) = P^n
π = πp 或πj = ∑(j=1~N)πiPij, j = 1,2,3, . ,N
有p(n) = p(0)P(n) -> p(0)P^n，且p(0) = π
-->得 πp^n=>πp^(n-1)=>πp=>π
即得出上述結論：鏈在任一時刻n∈T1的分布p(n)永遠與π一致
無意識統計學家法則
(Law Of The Unconscious Statistician, LOTUS)
在概率論和數理統計中，如果知道隨機變量X的概率分布，但是並不顯式地知道函數g(X)的分布，那麼LOTUS則是一個可以用來計算關於隨機變量X的函數g(X)的期望E[g(X)]的定理，該法則的具體形式依賴於隨機變量X概率分布的描述形式
如果隨機變量X的分布是離散的，而且我們知道他的PMF(概率質量函數)為fX，但不清楚fg(x),則無意識統計學家法則中計算g(X)的期望是：
E[g(X)] = ∑(x)g(x)fX(x)
其中的和式是在取遍X的所有可能值x之後求得的
如果隨機變量X的分布是連續的而且我們知道他的PDF(概率密度函數)為fX，但不知道fg(x)，那麼g(X)的期望為：
E[g(X)] = ∫(-∞->+∞)g(x)fX(x)
*總結：一直隨機變量X的概率分布，但是並不知道g(X)的分布，此時使用LOTUS公式法則能夠計算出函數g(X)的數學期望，其實就是在計算期望時，使用已知的X的PDF/PMF來代替未知的g(X)的g(X)的PDF/PMF
蒙特卡洛法求解定積分
1.投點法
投點法是蒙特卡洛法基本思想中最基礎也是最直觀的實例
如上圖所示，存在一個函數f(x)，若求解a->b的定積分(就是求解曲線下方的面積)
較容易算得面積的矩形罩在函數的積分區間上(上圖中的正方形，設其面積為Area)
然後隨機地向這個矩形框內投點，落在函數f(x)下方積分域內的為綠色點，其他的點為紅色點，之後統計下方積分域內的綠色點的數量佔全部點(紅點+綠點)數量的比例為r，之後就可以估算出函數f(x)從a->b的定積分為Area*r
hint：【蒙特卡洛法得出的值並不是精確值，而是近似值，該近似值隨著投點數量的增大而愈發精確 服從Week大數定理】
2.期望法(又稱平均值法)
任取一組相互獨立同分布的隨機變量{Xi}，Xi在[a, b]上服從分布律fX，也就是說，fX為隨機變量的PDF或PMF(概率質量函數/概率密度函數)
令g*(x) = g(x)/fX(x)，則g*(Xi)也是一組獨立同分布的隨機變量，而且因為g*(x)是關於x的函數，根據LOTUS可得
E[g*(Xi)] = ∫(a -> b)g*(x)fX(x)dx = ∫(a -> b)g(x)dx = I
由強大數定理
Pr(Lim(n->+∞)1/N∑(i=1~N)g*(Xi) = I) = 1
若選I^ = 1/N∑(i=1~N)g*(Xi)
則I^依據概率1收斂到I，平均值法使用I^作為I的近似值
假設我們要求的積分存在如下形式：
I = ∫(a -> b)g(x)dx
其中被積函數g(x)在區間[a, b]內可積，任意選擇一個簡便方法可以進行採樣的概率密度函數fX(x)，使其滿足下列條件：
(1).當g(x) != 0時，fX(x) ！= 0(a<=x<=b);
(2).∫(a -> b)fX(x)dx = 1;
|---->g(x)/fX(x), fX(x)!=0
記 g*(x) = |
| ---->0, fX(x) = 0
則原積分式可以寫為 I= ∫(a -> b)g*(x)fX(x)dx
因此求解積分的步驟為：
(1).產生服從分布律fX(x)的隨機變量Xi(i = 1,2,3, ., N)
(2).計算均值 I^ = 1/N∑(i=1~N)g*(Xi)
並將I^用作I的近似值，即I≈I^
如果a,b為有限值，則fX可取為均勻分布
fX(x) = 1/(b-a), a <= x <= b OR 0, otherwise
此時原來的積分式變為：
I = (b-a)∫(a -> b)g(x)*1/(b-a)dx
因此求解積分的步驟為：
(1).產生[a,b]上的均勻分布隨機變量Xi(i = 1,2,3, ., N)
(2).計算均值 I^ = (b-a)/N∑(i=1~N)g(Xi)
最後用I^作為I的近似值，即I≈I^
平均值法的幾何直觀解釋
注意積分的幾何意義為[a,b]區間內曲線下方的面積
當在[a,b]之間隨機取一點x時，對應的函數值就是f(x)，然後可以用f(x)*(b-a)粗略估算曲線下方的面積(即積分)
之後可以在區間[a,b]之間隨機取一系列點xi時(xi滿足均勻分布)，然後把估算出來的面積取平均作為積分估計得一個更好的近似值，因此得，採樣點越多則對於該區間的積分則越來越精確
按照上述思路，得到蒙特卡洛平均值法求解積分的最終積分公式為：
I^ = (b-a)1/N∑(i=1~N-1)f(Xi) = 1/N∑(i=1~N-1)f(Xi)/{1/(b-a)}
注意⚠️！：此處的1/(b-a)就是均勻分布的PMF(概率質量函數)
Hint 引入：定積分值誤差檢驗
方根誤差（Root Square Error）是觀測值與真值（理論值）偏差平方的平方根，是用來衡量觀測值同真值之間的偏差。隨採樣數量的增加，誤差逐漸降低。在達到一定數量的採樣之後，誤差變化不大。說明採樣次數提高到一定值後，對運算精度變化影響較少。進而採取改變抽樣法來減小方差，從而提高積分計算的精度
蒙特卡洛投點法應用
求圓周率pi的近似值：
(1)正方形內部有一個相切的圓，它們的面積之比是π/4。現在，在這個正方形內部，隨機產生10000-100000個點，計算它們與中心點的距離，從而判斷是否落在圓的內部,若這些點均勻分布，則圓周率 pi=res/n。
其中res:表示落到圓內投點數 n:表示總的投點數
(2)代碼[Python實現]
#!/usr/bin/env python # -*- coding:utf-8 -*-import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.patches import Circle# 投點次數10000-100000一萬次投擲到十萬次投擲# n = 100000n = 10000# 圓的半徑、圓心r = 1.0a,b = (0.,0.)# 正方形區域x_min, x_max = a-r, a+ry_min, y_max = b-r, b+r# 在正方形區域內隨機投點x = np.random.uniform(x_min, x_max, n) #均勻分布y = np.random.uniform(y_min, y_max, n)#計算點到圓心的距離d = np.sqrt((x-a)**2 + (y-b)**2)#統計落在圓內點的數目res = sum(np.where(d < r, 1, 0))#計算pi的近似值（Monte Carlo:用統計值去近似真實值）pi = 4 * res / nprint('pi: ',pi)#畫個圖fig = plt.figure()axes = fig.add_subplot(111)axes.plot(x, y,'ro',markersize = 1)plt.axis('equal')circle = Circle(xy=(a,b), radius=r ,alpha=0.5)axes.add_patch(circle)plt.show()[一萬次投擲模擬結果PI = 3.1556 圓的半徑為1]
[十萬次投擲模擬結果PI =3.15284 圓的半徑為10]
蒙特卡洛平均值法應用
積分函數為：y=3*x^2+4*cos(x)-4*x*sin(x)
import numpy as npimport pylab as pl
def fy(x): return x**3 + 4*x*np.cos(x) def dfy(x): return 3*x**2+4*np.cos(x)-4*x*np.sin(x)
def monto(x,a,b): return (b-a)/len(x)*sum(dfy(x))
def integral(a,b): return fy(b) - fy(a)
a,b,n=10,15,1000vm=np.ones(n)X=np.linspace(a,b,n)y=np.ones(n)vi=integral(a,b)for i in range(n): x = np.random.uniform(a,b,i+1) vm[i] = monto(x,a,b) y[i] = dfy(X[i])
pl.subplot(212)pl.title(" RSE ")pl.plot(range(n), np.sqrt((vi-vm)**2) )pl.xlabel("Monte Carlo sampling point")
pl.subplot(211)pl.title("function curve")pl.plot(X,y,label='$y=3*x^2+4*cos(x)-4*x*sin(x)$')pl.legend(loc='upper left')pl.show()蒙特卡洛採樣
在上述的例子中演示了蒙特卡洛發的基本思想，其中的核心思想就是使用隨機數，當所求解的問題可以轉化為某種隨機分布的特徵數(例如隨機事件出現的概率或者隨機變量的期望值等)時，往往可以考慮採用蒙特卡洛法。通過隨機採樣的方法，以隨機事件出現的頻率估計其概率，或者以採樣的數字特徵估算隨機變量的數字特徵，並將其作為問題的解，此種方法常用於求解複雜的高維積分問題
在實際的應用中，所要面對的第一個問題就是如何進行數據採樣，在計算機模擬中，通常所說的採樣指的是：從一個概率分布中生成觀察值的方法，且該分布通常是由它的概率密度函數決定並表示，但是結合之前所提到的無意識統計學家法則LOTUS，即使在已知PDF的情況下，讓計算機自動生成觀測值也並非易事，本質上來說，計算機只能模擬實現對均勻分布的採樣
(1).逆採樣方法(Inverse Sampling)
比較簡單的一種情況是，可以直接通過PDF(概率密度函數)與CDF(累積分布函數)之間的關係，求出對應的CDF，或者根本不知道PDF(概率密度函數)，但是知道CDF(累積分布函數)，此時就可以直接使用CDF的反函數(以及分位數函數)的方法進行採樣，此種方法又被稱為逆變換採樣方法(Inverse Transform Sampling)
逆採樣方法概述：假設已經得到了CDF的反函數F^-1(u)，如果想得到m個觀察值，則重複下列步驟m次即可：
(1).從最簡單常見易採樣的均勻分布Uniform(0, 1)中隨機生成一個值，並用u表示
(2).計算F^-1(u)的值x，則x就是從目標分布f(x)中的得出的一個有效採樣點
(3).重複上述過程m次
Hint：對於常見的均勻分布Uniform(0, 1)是非常容易採樣樣本的，一般通過線性同餘發生器可以很方便的生成(0, 1)之間的偽隨機數樣本。而其他常見的概率分布，無論是離散的分布還是連續的分布，它們的樣本都可以通過Uniform(0, 1)的樣本轉換而得
以指數分布為例，說明如何通過均勻分布來採樣服從指數分布的樣本集。指數分布的概率密度函數PDF為：
那麼它的概率分布函數CDF為：
下圖為指數分布和均勻分布的CDF圖，從左圖上看，在x ≥ 0的部分是一個單調遞增的函數(在定義域上單調非減)，定義域和值域是[0,+∞) -> [0,1)在p(x)大的地方它增長快，反之亦然
因為它是唯一映射的(在>0的部分，接下來我們只考慮這一部分)，所以它的反函數可以表示為F^-1(a)，a ∈ [0, 1)，值 域 為 [0, +∞)上圖的左圖就是F(x)的圖)因為F單調遞增，所以F^-1也是單調遞增的
利用反函數的定義，得到：
接下來定義[0, 1]均勻分布的CDF
根據公式(4), (5)得出：
因為F(x)的值域[0,1)根據公式(5), (6)可以改寫為：
根據F(x)的定義，它是exp分布的CDF，所以公式(7)的意思是F^-1(a)我們通過F的反函數將一個0到1均勻分布的隨機數轉換成了符合exp分布的隨機數，注意，以上推導對於CDF可逆的分布都是一樣的，對於exp來說，它的反函數的形式是：
[映射關係為：從y軸0-1的均勻分布樣本(green)映射得到服從指數分布的樣本(red)]
逆採樣方法代碼如下：
import matplotlib.pyplot as plt #python繪圖包def sampleExp(Lambda = 2,maxCnt = 50000): ys = [] standardXaxis = [] standardExp = [] for i in range(maxCnt): #循環產生50000個隨機數採樣 u = np.random.random() y = -1/Lambda*np.log(1-u) #F-1(X)反函數計算 ys.append(y) for i in range(1000): #循環 t = Lambda * np.exp(-Lambda*i/100) standardXaxis.append(i/100) standardExp.append(t) plt.plot(standardXaxis,standardExp,'r') plt.hist(ys,1000,normed=True) plt.show()
sampleExp()模擬採樣結果圖
隨機採樣樣本個數maxCnt=5000 Xaxis精度範圍0-5
隨機採樣樣本個數maxCnt=50000 Xaxis精度範圍0-10
隨機採樣樣本個數maxCnt=500000 Xaxis精度範圍0-10
隨機採樣樣本個數maxCnt=5000000 Xaxis精度範圍0-10
隨機採樣樣本個數maxCnt=5000000 Xaxis精度範圍0-100
由上述圖像可得：隨著採樣樣本數量的增加，概率直方圖越來越接近真實的CDF(累計分布函數)，而且採樣個數越多，與真實曲線越接近，曲線越光滑
拒絕採樣與自適應拒絕採樣
之前提到的逆變換採樣同樣也存在缺點，那就是有些分布的CDF可能難以通過對PDF的積分得到，再或者CDF的反函數根本不能輕易求出，為了彌補逆變換採樣的缺點，引申出了拒絕採樣(Reject Sampling)
設想對PDF為p(x)的函數進行採樣，但是由於種種原因(例如這個函數非常複雜)，對其進行採樣是相對困難的，但是另外一個PDF為q(x)的函數則相對的容易採樣，例如採用逆變換方法就可對其進行很容易的採樣，甚至q(x)本身就是一個均勻分布，那麼，當將q(x)與一個常數M相乘之後，可以實現下圖所示的關係，即M*q(x)將p(x)完全地「罩住」
然後重複下列步驟，知道m個被接受的採樣點：
(1). 從q(x)中獲得一個隨機採樣點xi;
(2). 對於xi計算接受概率(Acceptance Probability)
å = p(xi)/M*q(xi)
(3). 從Uniform(0, 1)中隨機生成一個值，u來表記
(4). 如果å>=u，則接受xi作為一個來自p(x)的採樣值，否則就拒絕xi並且回到步驟(1)重複過程
關於拒絕採樣的接受率問題
M的求取問題：很簡單，看圖有 Mq(x) >= p(x) {p(x)同π(x)}，那麼m>= p(x)/q(x)，求p(x)/q(x)的最大值 max[p(x)/q(x)]，即為m
簡要的從圖中直觀的解釋他的基本原理：在上圖中，哪些位置抽取出來的點較為容易被接受？很顯然，上方曲線Mq(x)和下方曲線π(x){同p(x)}所示的的間距最小，距離最接近的地方所被接受的概率較高，所以在這樣的範圍內所採樣的點會比較多，與此相反，兩個函數曲線間距較大的地方採樣的點會比較少，因為距離最遠的地方所被接受的概率最低，依據接受率，會被拒絕採樣 所以這樣就保證了拒絕採樣方法的有效性
拒絕採樣的應用實例
假設想對PDF：
0.3/np.sqrt(np.pi*2)*np.exp(-(x+2)**2/2) + 0.7/np.sqrt(np.pi*2)*np.exp(-(x-4)**2/2)的該函數進行拒絕採樣，具體實現代碼如下
#!/usr/bin/env python # -*- coding:utf-8 -*-import numpy as npimport matplotlib.pyplot as plt
def f1(x): # 假設f1為我們想要進行抽樣的分布 return 0.3/np.sqrt(np.pi*2)*np.exp(-(x+2)**2/2) + 0.7/np.sqrt(np.pi*2)*np.exp(-(x-4)**2/2) x = np.arange(-10, 10, 0.0001)plt.plot(x, f1(x), color='red', label='Target Density') # 畫出這個分布的圖像
# 計劃採用[-10, 10]之間的均勻分布作為建議分布size = int(1e+8)s = np.random.uniform(-10, 10, size) # 生成數據sy = 1/20k = 6 # 保證sy*k>f(x)plt.plot([-10, 10], [sy*k, sy*k], color='black') # 畫出建議分布的圖像t = 0.3/np.sqrt(np.pi*2)*np.exp(-(s+2)**2/2) + 0.7/np.sqrt(np.pi*2)*np.exp(-(s-4)**2/2)ran = np.random.uniform(0, k*sy, size)sampling = s[ran < t] # 生成的隨機數小於待抽樣分布的函數值則接受plt.hist(sampling, bins=125, normed=True, color='blue', label='Rejection Sampling')plt.legend(prop={'size': 20})plt.show()最終實現的拒絕採樣圖像為：
由上圖可以得知，圖中黑色線段為建議分布的K倍，在這裡指的是[−10, 10]之間的均勻分布的K倍，紅色線條表示待抽樣的分布，藍色部分是經過接受拒絕抽樣法得到的抽樣結果所服從的分布，該分布結果已經和原始待抽樣的分布大致契合，說明了接受拒絕抽樣法的有效性。
*拒絕採樣的弱點在於，當被拒絕的點很多時，採樣的效果很不理想，同時如果能夠找到一個跟目標分布函數非常接近的參考函數，那麼就可以保證被接受的點佔有大多數(被拒絕的點存在很少)
馬爾可夫鏈蒙特卡洛方法中的採樣方法
(Markov Chain Monte Carlo Sampling Method, MCMC Sampling Method)
回顧我們學過的概率論與數理統計中，在以貝葉斯定理為基礎的機器學習方法中，通常需要計算後驗概率，再通過最大化後驗概率(MAP)的方法進行參數推斷和決策
然而在很多時候，後驗分布的形式可能非常複雜，此時通過尋找最大化後驗概率估計或者對後驗概率進行積分等計算會十分困難，此時通過採樣方法中的，馬爾可夫蒙特卡洛法進行求解
重要性採樣(Importance Sampling)
重要性採樣中運用到了之前所提及到的蒙特卡洛平均值法 在計算定積分∫(-∞->+∞)g(x)dx時，會把g(x)拆解為兩項的乘積，即g(x) = f(x)p(x)，其中f(x)時某種形式的函數，而p(x)是關於隨機變量X的概率分布(PDF或者PMF)如此一來，上述定積分就可以表示為求f(x)的期望，即
∫(-∞->+∞)g(x)dx = ∫(-∞->+∞)f(x)p(x)dx = E[f(x)]
當g(x)的分布函數具有非常複雜的形式時，仍然無法直接求解，這時就可以用採樣的方式去近似，此時積分公式可以表示為
∫(-∞->+∞)g(x)dx = E[f(x)] = (1/n)∑(i=1~N)f(Xi)
在貝葉斯分析中，蒙特卡洛積分運算常常被用來在對後驗概率進行積分時做近似估計，例如要計算積分I(y) = ∫(a->b)f(y | x)p(x)dx時，可以使用近似形式I^(y) = (1/n)∑(i=1~N)f(y | xi)
在蒙特卡洛法進行積分求解時，非常重要的一個環節就是從特定的分布中進行採樣，這裡的「採樣」也即是生成某種滿足某種分布的觀測值，之前所言及的「逆採樣」和「拒絕採樣」方法同理
與之前的採樣方法(「逆採樣」/「拒絕採樣」)形成對比，「逆採樣，拒絕採樣」都是先對分布進行採樣，之後再用採樣的結果近似計算積分，而「重要性採樣(Importance Sampling)」為兩步並為一步，直接近似的計算求取積分
假設現在的目標為計算如下積分：
E[f(x)] = ∫f(x)p(x)dx
按照蒙塔卡羅球定積分的方法，將從滿足p(x)的概率分布中獨立地採樣出一系列隨機變量xi，然後便有E[f(x)] ≈ (1/n)∑(i=1~N)f(xi)
但是如果存在p(x)的形式非常複雜，對其的概率分布求解非常困難，此時可以做等量變換，有
∫f(x)p(x)dx = ∫f(x){p(x)/q(x)}q(x)dx
把其中的f(x){p(x)/q(x)}看作為一個新的函數h(x)則有
∫f(x)p(x)dx = ∫h(x)q(x)dx ≈ 1/n)∑(i=1~N)f(xi){p(xi)/q(xi)}
其中的{p(xi)/q(xi)}是xi的權重或稱為重要性權重(Importance Weight)，故稱該種採樣的方法為重要性採樣
[重要性採樣原理圖示]
如上圖所示，在使用重要性採樣時，並不會拒絕掉某些採樣點，與之前的「拒絕採樣不同」此時重要性採樣中的所有點均會使用，但是他們的權重是不同的，因為權重為{p(xi)/q(xi)}是xi的權重或稱為重要性權重(Importance Weight)，當p(xi) > q(xi)時，採樣點xi的權重大，反之則小，重要性採樣就是通過這樣一種方式從「參考分布」q(x)中獲取「目標分布」p(x)的方法
Hint 關於Importance Sampling(重要性採樣)的幾點看法及總結：關於重要性採樣，該方法與之前介紹的拒絕採樣法雖有不同，但實則換湯不換藥，方法名為「重要性採樣」實則採樣，但是最後得出的結論公式反而為求取均值，為何最後沒有得到樣本？
拒絕採樣時：通過接收判斷接受概率(Acceptance Probability)
å = p(xi)/M*q(xi)
來對通過q(x)得到的樣本進行篩選，使得最後得到的樣本分布律服從給定的p(x)，包絡的間隔越小接受程度越大，間隔越大接受程度越小
重要性採樣時：對於通過q(x)得到的樣本全部接受，不予拒絕，但是全部接受存在一個問題，之中得到的樣本點分布不回服從所給的p(x)分布，因此引入了重要性權重，來消除校正樣本全部接受帶來的偏差問題，參照上圖
p(z*)/q(z*)值越大的樣本，所給的權重越大，該值越小的樣本，所給的權重越小
最後，該重要性權重乘以樣本之後的結果即服從p(x)分布
關於重要性採樣的實例
首先定義函數f(x)，具體的函數定義及圖像如下：
# 定義f(x)def f_x(x): return 1/(1 + np.exp(-x))實例1 當p(x)與q(x)比較接近的情況下的重要性採樣實驗
具體代碼實現如下：
#!/usr/bin/env python # -*- coding:utf-8 -*-import numpy as npimport pandas as pdimport scipy.stats as statsimport seaborn as snsimport matplotlib.pyplot as plt"""Copyright (c)By JINZEWEI(金澤緯)."""
# 定義f(x)def f_x(x): return 1/(1 + np.exp(-x))
fig = plt.figure(figsize=(8,2))ax = fig.add_subplot(1,1,1)
# 繪製圖形x = np.linspace(0, 4, 100)ax.plot(x, f_x(x), 'b.', label=r'$\frac{1}{1+exp(-x)}$')
ax.legend(fontsize=17)#plt.show()
"""p(x)與q(x)比較接近，接下來定義分布p(x)和q(x)這兩個分布的均值非常接近，從每一個分布中抽樣3000個數據"""# q(x)mu_appro = 3sigma_appro = 1q_x = [np.random.normal(mu_appro, sigma_appro) for _ in range(3000)]# p(x)mu_target = 3.5sigma_target = 1p_x = [np.random.normal(mu_target, sigma_target) for _ in range(3000)]
fig = plt.figure(figsize=(10,6))ax = fig.add_subplot(1,1,1)
# 畫出兩個分布的圖像sns.distplot(p_x, label="distribution $p(x)$")sns.distplot(q_x, label="distribution $q(x)$")
plt.title("Distributions", size=16)plt.legend()
#計算E[f(x)],當x~p(x)時x_px = np.mean([f_x(i) for i in p_x])print(x_px)# 當x~q(x), E[f(x)] = mean(p(x)/q(x)*f(x))的計算結果p_pdf = stats.norm(mu_target, sigma_target)q_pdf = stats.norm(mu_appro, sigma_appro)
print(np.mean([p_pdf.pdf(i)/q_pdf.pdf(i) * f_x(i) for i in q_x]), np.var([p_pdf.pdf(i)/q_pdf.pdf(i) * f_x(i) for i in q_x]))plt.show()兩個接近的分布圖如下：
最終結果為：
實例2 當p(x)與q(x)相差較遠的情況下的重要性採樣實驗
具體的函數定義以及分布的代碼如下：
#!/usr/bin/env python # -*- coding:utf-8 -*-import numpy as npimport pandas as pdimport scipy.stats as statsimport seaborn as snsimport matplotlib.pyplot as plt
"""Copyright (c)By JINZEWEI(金澤緯)."""
# 定義f(x)def f_x(x): return 1/(1 + np.exp(-x))
# p(x)mu_target = 3.5sigma_target = 1p_x = [np.random.normal(mu_target, sigma_target) for _ in range(3000)]
# q(x)mu_appro = 1sigma_appro = 1q_x = [np.random.normal(mu_appro, sigma_appro) for _ in range(3000)]
fig = plt.figure(figsize=(10,6))ax = fig.add_subplot(1,1,1)
# 畫出兩個分布的圖像sns.distplot(p_x, color='red', label="distribution $p(x)$")sns.distplot(q_x, color='black', label="distribution $q(x)$")
plt.title("Distributions", size=16)plt.legend()
# 當x~p(x)時, E[f(x)]np.mean([f_x(i) for i in p_x])print(np.mean([f_x(i) for i in p_x]))
# 採樣數量=3000# 當x~q(x), E[f(x)] = mean(p(x)/q(x)*f(x))total = 3000p_x = [np.random.normal(mu_target, sigma_target) for _ in range(total)]q_x = [np.random.normal(mu_appro, sigma_appro) for _ in range(total)]
p_pdf = stats.norm(mu_target, sigma_target)q_pdf = stats.norm(mu_appro, sigma_appro)
print(np.mean([p_pdf.pdf(i)/q_pdf.pdf(i) * f_x(i) for i in q_x]), np.var([p_pdf.pdf(i)/q_pdf.pdf(i) * f_x(i) for i in q_x]))
# 採樣數量5000# 當x~q(x), E[f(x)] = mean(p(x)/q(x)*f(x))total = 5000p_x = [np.random.normal(mu_target, sigma_target) for _ in range(total)]q_x = [np.random.normal(mu_appro, sigma_appro) for _ in range(total)]
p_pdf = stats.norm(mu_target, sigma_target)q_pdf = stats.norm(mu_appro, sigma_appro)
print(np.mean([p_pdf.pdf(i)/q_pdf.pdf(i) * f_x(i) for i in q_x]), np.var([p_pdf.pdf(i)/q_pdf.pdf(i) * f_x(i) for i in q_x]))
# 當x~q(x), E[f(x)] = mean(p(x)/q(x)*f(x))# 採樣數量=10000total = 10000p_x = [np.random.normal(mu_target, sigma_target) for _ in range(total)]q_x = [np.random.normal(mu_appro, sigma_appro) for _ in range(total)]
p_pdf = stats.norm(mu_target, sigma_target)q_pdf = stats.norm(mu_appro, sigma_appro)
print(np.mean([p_pdf.pdf(i)/q_pdf.pdf(i) * f_x(i) for i in q_x]), np.var([p_pdf.pdf(i)/q_pdf.pdf(i) * f_x(i) for i in q_x]))
# 當x~q(x), E[f(x)] = mean(p(x)/q(x)*f(x))# 採樣數量=100000total = 100000p_x = [np.random.normal(mu_target, sigma_target) for _ in range(total)]q_x = [np.random.normal(mu_appro, sigma_appro) for _ in range(total)]
p_pdf = stats.norm(mu_target, sigma_target)q_pdf = stats.norm(mu_appro, sigma_appro)
print(np.mean([p_pdf.pdf(i)/q_pdf.pdf(i) * f_x(i) for i in q_x]), np.var([p_pdf.pdf(i)/q_pdf.pdf(i) * f_x(i) for i in q_x]))
plt.show()兩個距離較遠的分布圖如下：
最終結果為：
總結：由上述採樣容量以及方差的大小程度可知：兩個分布圖像差的比較多的時候，需要更多採樣來穩定mean(均值)較準確的估計，但與此同時估計的方差偏移比較大，因為來自於bias = p(x)/q(x)，兩個分布的差值比較大所以bias的值會比較大，最終導致方差較大
馬爾可夫鏈蒙特卡洛法的基本概念梳理
馬爾可夫蒙特卡洛方法是一類用於從一個概率分布中進行採樣的算法，該法以馬爾可夫鏈為基礎，而被構造的馬爾可夫鏈(Markov Chain)分布應該同需要的分布等價
與之相關的馬爾可夫原理 & 馬爾可夫鏈原理 & Chapman-Kolmogorov(查普曼-科爾莫戈羅夫)方程 & 轉移矩陣等基礎知識不做過多贅述，參考本文初始部分
⚠️再次注意：細緻平衡原理⚠️：
一條馬爾可夫鏈擁有平穩分布的一個充分必要條件是，對於任意兩個狀態i, j其滿足細緻平衡(Detailed Balance)<==> p(j->i)πi = p(i->j)πj
可以見得，此時(πP)j = ∑(i=1~N)πip(i->j) = ∑(i=1~N)πjp(j->i) = πj∑(i=1~N)p(j->i) = πj 故，π = πP，其明顯滿足平穩分布的條件，如果一條馬爾可夫鏈滿足細緻平衡，則稱他是具有可逆性的(可逆的)
[細緻平穩原理圖例資料參考]
[馬爾可夫鏈狀態轉移圖例]
最後總結關於MCMC(Markov Chain Monte Carlo)的基本思想：在拒絕採樣和重要性採樣中，當前生成的樣本點與之前生成的樣本點之間不存在直接關係(無聯繫性)，它們的採樣都是獨立進行的，然而MCMC是基於Markov Chain進行的採樣，結合之前提及到的馬爾可夫性，表明當前的樣本點生成與上一時刻的樣本點是有關的 結合下圖
假設當前時刻生成的樣本點為x，下一個時刻採樣到的x'的轉移概率為p(x|x')，或簡潔表記為p(x->x') 依據馬爾可夫鏈的遍歷性原理，最終的概率分布收斂至π(x)「Hint：鏈在任一時刻n∈T1的分布p(n)永遠與π一致」 我們希望這個過程如此持續進行下去，也就最終收斂得到要採樣的分布
MCMC採樣算法圖示：
MCMC方法代表實例算法
1.米特羅波利斯-黑斯廷斯算法(Metropolis-Hastings)
對於給定的概率分布，如果直接從中採樣比較困難，那麼M-H算法就是直接從中獲取一系列隨機樣本的MCMC方法，M-H算法的執行步驟如下
米特羅波利斯-黑斯廷斯算法的執行步驟為：先隨機指定一個初始的樣本點x0，u時來自均勻分布的一個隨機數，x*是來自p分布的樣本點，該樣本點遵循MCMC方法是否從前一個位置跳轉到當前的x*由å(x*)和u的大小關係決定，如果接受則跳轉，即新的採樣點為x*，否則就拒絕，即新的採樣點仍然為上一輪的採樣點
證明M-H算法是否滿足細緻平衡有如下證明過程：
吉布斯採樣
吉布斯採樣可以看作是M-H算法的一個特例，該特殊之處在於，使用吉布斯採樣時，通常是對多變量進行採樣
吉布斯採樣的算法流程描述如下：
米特羅波利斯-黑斯廷斯算法與吉布斯採樣算法的實現
我們之前提到狀態轉移函數T的選擇，可以看到如果我們選擇一個對稱的轉移函數，即T(a,b)=T(b,a) 上面的接受概率還可以簡化為
這也是一般Metropolis算法中採用的方法，T使用一個均勻分布即可，所有狀態之間的轉移概率都相同。實驗中我們使用二元高斯分布等均勻分布來進行採樣模擬
代碼示例如下：
#!/usr/bin/env python # -*- coding:utf-8 -*-from matplotlib import pyplot as pltimport numpy as npfrom mpl_toolkits.mplot3d import Axes3Dfrom matplotlib import cmimport math#PI = 3.14159265358979323846PI = math.pidef domain_random(): return np.random.random()*3.8-1.9def get_p(x): #return 1/(2*PI)*np.exp(- x[0]**2 - x[1]**2) return np.sin(x[0]**2 + x[1]**2) + 1 #return 1/(2*PI*math.sqrt(1-0.25))*np.exp(-1/(2*(1-0.25))*(x[0]**2 -x[0]*x[1] + x[1]**2))
def get_tilde_p(x): return get_p(x)*20
def partialSampler(x,dim): xes = [] for t in range(10): xes.append(domain_random()) tilde_ps = [] for t in range(10): tmpx = x[:] tmpx[dim] = xes[t] tilde_ps.append(get_tilde_p(tmpx))
norm_tilde_ps = np.asarray(tilde_ps)/sum(tilde_ps) u = np.random.random() sums = 0.0 for t in range(10): sums += norm_tilde_ps[t] if sums>=u: return xes[t]
def plotContour(plot = False): X = np.arange(-2, 2, 0.05) Y = np.arange(-2, 2, 0.05) X, Y = np.meshgrid(X, Y) Z = get_p([X,Y]) plt.contourf(X, Y, Z, 100, alpha = 1.0, cmap =cm.coolwarm) plt.contour(X, Y, Z, 10, colors = 'purple', linewidth = 0.01)
if plot: plt.show()
def plot3D(): X = np.arange(-2, 2, 0.05) Y = np.arange(-2, 2, 0.05) X, Y = np.meshgrid(X, Y) Z = get_p([X,Y]) fig = plt.figure() ax = Axes3D(fig) ax.grid(False) ax.w_yaxis.set_pane_color((1,1,1,0)) ax.w_xaxis.set_pane_color((1,1,1,1)) ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm) plt.show()
plotContour()plot3D()
def metropolis(x): new_x = (domain_random(),domain_random()) acc = min(1,get_tilde_p((new_x[0],new_x[1]))/get_tilde_p((x[0],x[1]))) u = np.random.random() if u<acc: return new_x return x
def testMetropolis(counts ,drawPath = False): plotContour() x = (domain_random(),domain_random()) xs = [x] name = []
for i in range(counts): xs.append(x) x = metropolis(x) name.append(i)
if drawPath: plt.plot(list(map(lambda x:x[0],xs)), list(map(lambda x:x[1], xs)), 'k-', linewidth=0.05)
plt.scatter(list(map(lambda x:x[0], xs)), list(map(lambda x:x[1], xs)), c='black', s=10, marker='*') """for i in range(counts): plt.annotate(name[i], xy = (list(map(lambda x:x[0], xs))[i], list(map(lambda x:x[1], xs))[i]), xytext=(list(map(lambda x:x[0], xs))[i] + 0.1, list(map(lambda x:x[1], xs))[i] + 0.1), fontsize=5) """ plt.show() pass
def gibbs(x): rst = np.asarray(x)[:] path = [(x[0],x[1])] for dim in range(2): new_value = partialSampler(rst,dim) rst[dim] = new_value path.append([rst[0],rst[1]]) return rst,path
def testGibbs(counts = 100,drawPath = False): plotContour()
x = (domain_random(),domain_random()) xs = [x] paths = [x] name = [] for i in range(counts): xs.append([x[0],x[1]]) x,path = gibbs(x) paths.extend(path) name.append(i) if drawPath: plt.plot(list(map(lambda x:x[0], paths)), list(map(lambda x:x[1],paths)), c='k-', linewidth=0.05) plt.scatter(list(map(lambda x:x[0], xs)), list(map(lambda x:x[1], xs)), c='blue', s=5, marker='*') """for i in range(counts): plt.annotate(name[i], xy = (list(map(lambda x:x[0], xs))[i], list(map(lambda x:x[1], xs))[i]), xytext=(list(map(lambda x:x[0], xs))[i] + 0.1, list(map(lambda x:x[1], xs))[i] + 0.1), fontsize=5) """ plt.show() pass
testMetropolis(10, False)testMetropolis(100, False)testMetropolis(1000, False)testMetropolis(5000, False)
testGibbs(10, False)testGibbs(100, False)testGibbs(1000, False)testGibbs(5000, False)使用米特羅波利斯-黑斯廷斯算法&吉布斯採樣進行模擬採樣二元高斯分布的情形(Python代碼執行結果)
二元高斯函數圖像
[二元高斯函數圖像 平面等高線形式]
[使用M-H算法採樣二元高斯函數結果，樣本容量n=10]
[使用M-H算法採樣二元高斯函數結果，樣本容量n=100]
[使用M-H算法採樣二元高斯函數結果，樣本容量n=1000]
[使用M-H算法採樣二元高斯函數結果，樣本容量n=5000]
[使用吉布斯採樣算法採樣二元高斯函數結果，樣本容量n=10]
[使用吉布斯採樣算法採樣二元高斯函數結果，樣本容量n=100]
[使用吉布斯採樣算法採樣二元高斯函數結果，樣本容量n=1000]
[使用吉布斯採樣算法採樣二元高斯函數結果，樣本容量n=5000]
使用米特羅波利斯-黑斯廷斯算法&吉布斯採樣進行模擬採樣
sin[(x^2+y^2)^1/2]分布的情形(Python代碼執行結果)
sin[(x^2+y^2)^1/2]函數圖像
[sin[(x^2+y^2)^1/2]函數圖像 平面等高線形式]
[使用M-H算法採樣sin[(x^2+y^2)^1/2]函數結果，樣本容量n=10]
[使用M-H算法採樣sin[(x^2+y^2)^1/2]函數結果，樣本容量n=100]
[使用M-H算法採樣sin[(x^2+y^2)^1/2]函數結果，樣本容量n=1000]
[使用M-H算法採樣sin[(x^2+y^2)^1/2]函數結果，樣本容量n=5000]
[使用吉布斯採樣算法採樣sin[(x^2+y^2)^1/2]函數結果，樣本容量n=10]
[使用吉布斯採樣算法採樣sin[(x^2+y^2)^1/2]函數結果，樣本容量n=100]
[使用吉布斯採樣算法採樣sin[(x^2+y^2)^1/2]函數結果，樣本容量n=1000]
[使用吉布斯採樣算法採樣sin[(x^2+y^2)^1/2]函數結果，樣本容量n=5000]
最後
我使用了樣本容量為250的M-H採樣算法和吉布斯採樣算法來對比兩者之間的異同(參照採樣順序序號)
[M-H Sampling Size 250 採樣結果圖例]
[Gibbs Sampling Size 250 採樣結果圖例]
從採樣路徑來看，依照採樣順序序號：
*Metropolis是完全隨機的，以一個概率進行拒絕
*而Gibbs Sampling則是在某個維度上進行轉移(某連續時刻段中相當於在等高線相等的位置連續採樣等高位置不變)
整理不易，且學且珍惜
數學是一種精神，一種理性的精神。正是這種精神，激發、促進、鼓舞並驅使人類的思維得以運用到最完善的程度，亦正是這種精神，試圖決定性地影響人類的物質、道德和社會生活；試圖回答有關人類自身存在提出的問題；努力去理解和控制自然；盡力去探求和確立已經獲得知識的最深刻的和最完美的內涵。——克萊因《西方文化中的數學》
參考文獻：概率論與數理統計(浙江大學第四版)
機器學習中的數學修煉 左飛(清華大學出版社)
AI技術交流 + 興趣討論：