UA MATH575B 數(shù)值分析下 計(jì)算統(tǒng)計(jì)物理例題2

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

一個(gè)位于原點(diǎn)$x_0=0$的粒子源，它釋放出的例子下一個(gè)時(shí)刻有2/3的概率的停在原地，有1/3的概率往右移動(dòng)一個(gè)單位；當(dāng)粒子位于$x=1,2,?,19$的位置時(shí)，有2/3的概率往左移動(dòng)一個(gè)單位，有1/3的概率往右移動(dòng)一個(gè)單位；$x=20$是一個(gè)吸收狀態(tài)。求這個(gè)吸收狀態(tài)的吸收率$\gamma$。

理論解法

這個(gè)粒子的運(yùn)動(dòng)可以用一個(gè)Markov鏈表示，記概率轉(zhuǎn)移矩陣為$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)$表示粒子從狀態(tài)$x_0$出發(fā)轉(zhuǎn)移到狀態(tài)子集$A$中所需的時(shí)間的期望，則
$$E(x\to A)=\sum _{y\in A}T(x,y)+\sum _{y\in /A}T(x,y)(1+E(y\to A))$$
第一項(xiàng)表示經(jīng)過一步轉(zhuǎn)移就到了狀態(tài)子集$A$中的概率，第二項(xiàng)表示如果第一步?jīng)]能直接轉(zhuǎn)移到位，記轉(zhuǎn)移到了狀態(tài)$y$，對(duì)應(yīng)的概率為$T(x,y),y\in /A$，然后$E(y\to A)$表示粒子從狀態(tài)$y$出發(fā)轉(zhuǎn)移到狀態(tài)子集$A$中所需的時(shí)間的期望，但此時(shí)已經(jīng)走了一步了所以加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$$
當(dāng)$x=1,?,18$時(shí)，
$$E_x=1+\frac{2}{3}{E}_{x?1}+\frac{1}{3}{E}_{x+1}$$
當(dāng)$x=19$時(shí)，
$$E_{19}=1+\frac{2}{3}{E}_{18}$$
由此我們得到了一個(gè)有二十個(gè)未知量、二十個(gè)方程的線性方程組，求解這個(gè)線性方程組我們可以得到$E_0=6291390$。也就是說在原點(diǎn)放出一個(gè)粒子，它平均需要移動(dòng)6291390步才能到達(dá)吸收狀態(tài)，因此吸收狀態(tài)的吸收率是$\gamma =1/6291390$。

特征值法（近似解）

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

模擬解法

Rejection Sampling

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

根據(jù)這個(gè)code可以畫出停時(shí)的分布

估計(jì)停時(shí)的均值為6.2546e+6，吸收率是這個(gè)均值的倒數(shù)。這個(gè)采樣最大的問題是慢，因?yàn)橥易叩母怕矢　?/p>

Importance Sampling

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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。