如何創建pymc3容器

我試圖建立（在調和勢布朗粒子）的朗之萬方程的特定結果的似然函數的模型：如何創建pymc3容器

這裏是我的模型pymc2似乎工作： https://github.com/hstrey/BayesianAnalysis/blob/master/Langevin%20simulation.ipynb

#define the model/function to be fitted. 
def model(x): 
    t = pm.Uniform('t', 0.1, 20, value=2.0) 
    A = pm.Uniform('A', 0.1, 10, value=1.0) 

    @pm.deterministic(plot=False) 
    def S(t=t): 
     return 1-np.exp(-4*delta_t/t) 

    @pm.deterministic(plot=False) 
    def s(t=t): 
     return np.exp(-2*delta_t/t) 

    path = np.empty(N, dtype=object) 

    path[0]=pm.Normal('path_0',mu=0, tau=1/A, value=x[0], observed=True) 
    for i in range(1,N): 
     path[i] = pm.Normal('path_%i' % i, 
         mu=path[i-1]*s, 
         tau=1/A/S, 
         value=x[i], 
         observed=True) 
     return locals() 

mcmc = pm.MCMC(model(x)) 
mcmc.sample(20000, 2000, 10)

的基本思想是每個點取決於鏈（Markov鏈）之前的點。順便說一句，x是一個數組數組，N是它的長度，delta_t是時間步長= 0.01。任何想法如何在pymc3中實現這個？我試過了：

# define the model/function for diffusion in a harmonic potential 
DHP_model = pm.Model() 
with DHP_model: 
    t = pm.Uniform('t', 0.1, 20) 
    A = pm.Uniform('A', 0.1, 10) 

    S=1-pm.exp(-4*delta_t/t) 

    s=pm.exp(-2*delta_t/t) 

    path = np.empty(N, dtype=object) 

    path[0]=pm.Normal('path_0',mu=0, tau=1/A, observed=x[0]) 
    for i in range(1,N): 
     path[i] = pm.Normal('path_%i' % i, 
         mu=path[i-1]*s, 
         tau=1/A/S, 
         observed=x[i])

不幸的是，只要我嘗試運行它，模型就崩潰了。我在我的機器上嘗試了一些pymc3示例（教程），這是行得通的。

在此先感謝。我真的希望pymc3中的新採樣器能夠幫助我實現這個模型。我試圖將貝葉斯方法應用於單分子實驗。

來源

2016-04-15 Helmut Strey

您可以通過自定義分佈（通過擴展Continuous）來知道計算整個路徑的對數似然的公式，而不是在循環中創建許多單獨的正態分佈的一維變量。你可以從pymc3已經知道的正態似然公式中引導這個似然公式。一個例子見the built-in AR1 class。

因爲你的粒子遵循馬爾科夫特性，你可能看起來像

import theano.tensor as T 

def logp(path): 
    now = path[1:] 
    prev = path[:-1] 

    loglik_first = pm.Normal.dist(mu=0., tau=1./A).logp(path[0]) 
    loglik_rest = T.sum(pm.Normal.dist(mu=prev*ss, tau=1./A/S).logp(now)) 
    loglik_final = loglik_first + loglik_rest 

    return loglik_final

我猜，你想在每一個時間步長繪製ss的值，在這種情況下，你應該確保指定ss = pm.exp(..., shape=len(x)-1)，以便上面的塊中的prev*ss被解釋爲元素方式的乘法。

然後，你可以用

path = MyLangevin('path', ..., observed=x)

指定你的觀察這應該運行多更快。

來源

2016-08-31 04:35:39 Noah

由於我沒有看到我的問題的答案，讓我自己來回答。我提出了以下解決方案：

# now lets model this data using pymc 
# define the model/function for diffusion in a harmonic potential 
DHP_model = pm.Model() 
with DHP_model: 
    D = pm.Gamma('D',mu=mu_D,sd=sd_D) 
    A = pm.Gamma('A',mu=mu_A,sd=sd_A) 

    S=1.0-pm.exp(-2.0*delta_t*D/A) 

    ss=pm.exp(-delta_t*D/A) 

    path=pm.Normal('path_0',mu=0.0, tau=1/A, observed=x[0]) 
    for i in range(1,N): 
     path = pm.Normal('path_%i' % i, 
         mu=path*ss, 
         tau=1.0/A/S, 
         observed=x[i]) 

    start = pm.find_MAP() 
    print(start) 
    trace = pm.sample(100000, start=start)

不幸的是，這段代碼需要N = 50，在6小時到2天之間的任何地方編譯。我運行在運行Ubuntu的相當快的PC（24Gb RAM）上。我試圖使用GPU，但運行速度稍慢。我懷疑內存問題，因爲它在運行時使用了99.8％的內存。我嘗試了與斯坦同樣的計算，只需要2分鐘即可運行。

來源

2016-07-27 22:10:50

如何創建pymc3容器

回答

相關問題