UA MATH575B 數值分析下 計算統計物理例題2

• 理論解法
• • C-K方程法
  • 特征值法（近似解）
• 模擬解法
• • Rejection Sampling
  • Importance Sampling

一個位于原點$x_0=0$的粒子源，它釋放出的例子下一個時刻有2/3的概率的停在原地，有1/3的概率往右移動一個單位；當粒子位于$x=1,2,?,19$的位置時，有2/3的概率往左移動一個單位，有1/3的概率往右移動一個單位；$x=20$是一個吸收狀態。求這個吸收狀態的吸收率$\gamma$。

理論解法

這個粒子的運動可以用一個Markov鏈表示，記概率轉移矩陣為$T$：

%% Construct transition matrix T = zeros(21,21); T(1,1) = 2/3; T(21,21) = 1; for i = 1:21for j = 1:21if j == (i-1) && i<21T(i,j) = 2/3;elseif j == (i+1)T(i,j) = 1/3;endend end

C-K方程法

用Chapman-Kolmogorov方程的思想，記$E(x_0\to A)$表示粒子從狀態$x_0$出發轉移到狀態子集$A$中所需的時間的期望，則
$$E(x\to A)=\sum _{y\in A}T(x,y)+\sum _{y\in /A}T(x,y)(1+E(y\to A))$$
第一項表示經過一步轉移就到了狀態子集$A$中的概率，第二項表示如果第一步沒能直接轉移到位，記轉移到了狀態$y$，對應的概率為$T(x,y),y\in /A$，然后$E(y\to A)$表示粒子從狀態$y$出發轉移到狀態子集$A$中所需的時間的期望，但此時已經走了一步了所以加1，
$$\sum _{y\in A}T(x,y)+\sum _{y\in /A}T(x,y)(1+E(y\to A))=1+\sum _{y\in /A}T(x,y)E(y\to A)$$
記$E_x=E(x\to \{20\})$，則
$$E_0=1+\frac{2}{3}{E}_0+\frac{1}{3}{E}_1$$
當$x=1,?,18$時，
$$E_x=1+\frac{2}{3}{E}_{x?1}+\frac{1}{3}{E}_{x+1}$$
當$x=19$時，
$$E_{19}=1+\frac{2}{3}{E}_{18}$$
由此我們得到了一個有二十個未知量、二十個方程的線性方程組，求解這個線性方程組我們可以得到$E_0=6291390$。也就是說在原點放出一個粒子，它平均需要移動6291390步才能到達吸收狀態，因此吸收狀態的吸收率是$\gamma =1/6291390$。

特征值法（近似解）

假設初始狀態的概率分布是${\pi }_0$，則$n$步后的狀態分布是${\pi }_n={\pi }_0{T}^n$。做轉移概率矩陣的對角化：$\hat{T}=V\Lambda V^{?1}$，則${\pi }_n={\pi }_{0}V{\Lambda }^n{V}^{?1}$。注意到這個Markov鏈只有一個吸收狀態，所以不論${\pi }_0$取什么值，當$n$足夠大時，${\pi }_n$都會退化為$\delta (x_n?20)$，并且狀態$20$對應的特征值為0。因此我們可以用第二大的特征值來近似${\pi }_n$，
$${\pi }_n\approx {\pi }_0{\lambda }_2^n$$
吸收率就近似為$1?{\lambda }_2$

模擬解法

Rejection Sampling

這是老師給的答案，大概是用python來寫，先import一下random的包然后設置初始值。當粒子沒有到狀態20時做rejection sampling，輸出所有例子最快到達吸收狀態的時間。

根據這個code可以畫出停時的分布

估計停時的均值為6.2546e+6，吸收率是這個均值的倒數。這個采樣最大的問題是慢，因為往右走的概率更小。

Importance Sampling

在rejection sampling中，我們用的是原始的Markov鏈的轉移概率，現在考慮做重要性采樣加快模擬的速度。這個Markov鏈比較有意思，相當于是一個不對稱的一維的隨機游走，所以要構造一個新的Markov鏈來采樣只需要一個參數$p$（code第二行，手動輸入）。接下來用這個新的Markov鏈做rejection sampling，$x$記錄的是新的Markov鏈的狀態轉移，compensation記錄的是importance sampling的權重，因此每次計數需要加compensation乘以1。
從這個結果可以看出，$p$越大吸收率估計量扭曲得越厲害（Markov鏈被扭曲了），但速度會加快。要在扭曲和速度之間做tradeoff，可以把這兩個Markov鏈合起來，當狀態小于5時就用原來的Markov鏈，大于5時就用參數為$p$的Markov鏈。

總結

以上是生活随笔為你收集整理的UA MATH575B 数值分析下 计算统计物理例题2的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。