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

• 統計物理方法的解析解
• Markov鏈
• • 理論解
  • 數值解
• Monte Carlo模擬.

一道有趣的統計物理的題目。下面這個簡單的迷宮中，一只老鼠一開始在3的位置，它在岔路口走每條路的概率是均等的，比如這個時刻在3的位置，那么下一個時刻有1/2的可能在2的位置，有1/2的時刻在7的位置，目標是計算老鼠成功逃出的概率。

統計物理方法的解析解

首先我們可以試圖找一下理論解。不妨假設3的位置是一塊榴蓮，13的位置是一個效果非常好的活性炭除臭盒，12的位置會讓氣味集中，我們想知道榴蓮的味道均勻擴散后（達到平衡狀態）迷宮里面臭味的分布情況。根據平衡狀態的關系，假設用$r_i$表示平衡狀態每個位置的臭味，以位置0為例，
$$r_0=(r_1+r_4)/2$$
除位置12、13外一共有14個這樣的方程（因為12為1，13為0），可以用mathematica來解（默默吐槽一下，這個就是算平穩分布而已）：

這些值，比如$r_3=9/23$，說明位置3的氣味分子會有$9/23$流竄到位置12。也就是說小鼠從位置3出發逃出去的概率是9/23=0.391304348…。

Markov鏈

寫出上面那個迷宮的Markov轉移概率矩陣：

理論解

假設初始狀態的概率分布是${\pi }_0$，則$n$步后的狀態分布是${\pi }_n={\pi }_0{\hat{T}}^n$。做轉移概率矩陣的對角化：$\hat{T}=V\Lambda V^{?1}$，則${\pi }_n={\pi }_{0}V{\Lambda }^n{V}^{?1}$。注意到這個Markov鏈有兩個吸收狀態（位置12和13），假設平穩分布是$\pi$，則$\pi =\pi \hat{T}$，也就是說吸收狀態對應的轉移概率矩陣的特征值為1，也就是說$\lim {\Lambda }^n(13)=1,\lim {\Lambda }^n(14)=1$。所以
$${\pi }_{\infty }(13)=(V?V^{?1})(13),{\pi }_{\infty }(14)=(V?V^{?1})(14)$$
用python求解

可以得到概率大約是是0.39130435。

數值解

數值計算的思路是用${\pi }_n={\pi }_0{\hat{T}}^n$做矩陣乘法，當$n$足夠大的時候${\pi }_n$會收斂到平穩分布。


這里估計的概率是0.391304348。基于Markov鏈的這兩種方法都還是很準確的。

Monte Carlo模擬.

一個做Monte Carlo模擬的matlab示例。模擬百萬次成功逃出的次數是390878，所以估計概率是0.390878。

P = [0 0.5 0 0 0.5 0 0 0 0 0 0 0 0 0 0 0;1/3 0 1/3 0 0 1/3 0 0 0 0 0 0 0 0 0 0; ...0 0.5 0 0.5 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0.5 0 0 0 0 0.5 0 0 0 0 0 0 0 0; ...1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0.5 0 0 0 0 0 0 0 0.5 0 0 0 0 0 0; ...0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0; 0 0 0 1/3 0 0 1/3 0 0 0 0 1/3 0 0 0 0; ...0 0 0 0 0 0 0 0 0 0.5 0 0 0.5 0 0 0; 0 0 0 0 0 1/3 0 0 1/3 0 1/3 0 0 0 0 0; ...0 0 0 0 0 0 0 0 0 0.5 0 0 0 0 0.5 0; 0 0 0 0 0 0 0 0.5 0 0 0 0 0 0 0 0.5; ...0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0; ...0 0 0 0 0 0 0 0 0 0 1/3 0 0 1/3 0 1/3; 0 0 0 0 0 0 0 0 0 0 0 0.5 0 0 0.5 0];rng(0); Eaten = 0; Escape = 0; for i = 1:1000000 % Set free 1 million micex = 3; % Set free 1 million mice at position 3while x ~= 12 && x ~= 13 % when the mouse hasn't escaped or been eaten u = rand(1,1);sum = P(x+1,1);y = 0;while sum < u % check which state will y transit tosum = sum + P(x+1,y+2);y = y + 1;endx = y;endswitch xcase 12Eaten = Eaten + 1; continue; % record the number of eaten micecase 13Escape = Escape + 1; continue; % record the number of escaped miceend end

總結

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

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