在pymc3中重寫動態系統中參數估計的pymc腳本

Question

我想使用pymc3來估計Hodgkin Huxley神經元模型中的未知參數和狀態。我在pymc中的代碼基於http://healthyalgorithms.com/2010/10/19/mcmc-in-python-how-to-stick-a-statistical-model-on-a-system-dynamics-model-in-pymc/，執行得相當好。

#parameter priors 
@deterministic 
def HH(priors in here) 
    #model equations 
    #return numpy arrays that somehow contain the probability distributions as elements 
    return V,n,m,h 

#Make V deterministic in one line. Seems to be the magic that makes this work. 
V = Lambda('V', lambda HH=HH: HH[0]) 

#set up the likelihood 
A = Normal('A',mu=V,tau=2.0,value=voltage_data,observed = True) 
#run the sampling... 

在pymc3中，Lambda技巧並不適用於我。這裏是我的嘗試之一：

import numpy as np 
import theano.tensor as tt 
from pymc3 import Model, Normal, Uniform, Deterministic, sample, NUTS, Metropolis, find_MAP 
import scipy 

#observed data 
T = 10 
dt = 0.02 

voltage_data_file = 'noise_measured.dat' 
voltage_data = np.loadtxt(voltage_data_file) 
voltage_data = voltage_data[0:T] 

current_data_file = 'current.dat' 
current_data = np.loadtxt(current_data_file) 
current_data = current_data[0:T] 

#functions needed later 
def x_inf(V,vx,dvx): 
    return 0.5*(1 + np.tanh((V - vx)/dvx)) 
def tau(V,vx_t,dvx_t,tx_0,tx_1): 
    return tx_0 + 0.5*tx_1*(1 + np.tanh((V- vx_t)/dvx_t)) 

#Model 
NaKL_Model = Model() 
with NaKL_Model: 
    #priors 
    g_Naa = Uniform('g_Naa',0.0,150.0) 
    E_Na = Uniform('E_Na',30.0,70.0) 
    g_K = Uniform('g_K',0.0,150.0) 
    E_K = Uniform('E_K',-100.0,-50.0) 
    g_L = Uniform('g_L',0.1,1.0) 
    E_L = Uniform('E_L',-90.0,-50.0) 
    vm = Uniform('vm',-60.0,-30.0) 
    vm_t = Uniform('vm_t',-60.0,-30.0) 
    dvm = Uniform('dvm',10.0,20.0) 
    dvm_t = Uniform('dvm_t',10.0,20.0) 
    tmm_0 = Uniform('tmm_0',0.05,0.25) 
    tm_1 = Uniform('tm_1',0.1,1.0) 
    vh = Uniform('vh',-80.0,-40.0) 
    vh_t = Uniform('vh_t',-80.0,-40.0) 
    dvh = Uniform('dvh',-30.0,-10.0) 
    dvh_t = Uniform('dvh_t',-30.0,-10.0) 
    th_0 = Uniform('th_0',0.5,1.5) 
    th_1 = Uniform('th_1',5.0,9.0) 
    vn = Uniform('vn',-70.0,-40.0) 
    vn_t = Uniform('vn_t',-70.0,-40.0) 
    dvn = Uniform('dvn',10.0,40.0) 
    dvn_t = Uniform('dvn_t',10.0,40.0) 
    tn_0 = Uniform('tn_0',0.5,1.5) 
    tn_1 = Uniform('tn_1',3.0,7.0) 
    #Initial Conditions 
    V_0 = Uniform('V_0',-100.0,50.0) 
    n_0 = Uniform('n_0',0.0,1.0) 
    m_0 = Uniform('m_0',0.0,1.0) 
    h_0 = Uniform('h_0',0.0,1.0) 

    def HH(i,V_current,n_current,m_current,h_current,g_Naa=g_Naa,E_Na=E_Na,g_K=g_K,E_K=E_K,g_L=g_L,E_L=E_L,vm=vm,vm_t=vm_t,dvm=dvm,dvm_t=dvm_t,tmm_0=tmm_0,tm_1=tm_1,vh=vh,vh_t=vh_t,dvh=dvh,dvh_t=dvh_t,th_0=th_0,th_1=th_1,vn=vn,vn_t=vn_t,dvn=dvn,dvn_t=dvn_t,tn_0=tn_0,tn_1=tn_1): 

     V_new = V_current + dt*(g_L*(E_L - V_current) + g_K*n_current**4*(E_K - V_current) + g_Naa*m_current**3*h_current*(E_Na - V_current) + current_data[i]) 

     n_new = n_current + dt*(x_inf(V_current,vn,dvn)-n_current)/tau(V_current,vn_t,dvn_t,tn_0,tn_1) 
     m_new = m_current + dt*(x_inf(V_current,vm,dvm)-m_current)/tau(V_current,vm_t,dvm_t,tmm_0,tm_1) 
     h_new = h_current + dt*(x_inf(V_current,vh,dvh)-h_current)/tau(V_current,vh_t,dvh_t,th_0,th_1) 

     return V_new,n_new,m_new,h_new 

    V_state = [V_0] 
    n_state = [n_0] 
    m_state = [m_0] 
    h_state = [h_0] 

    #A = [Normal('A0',mu=V_0,tau=2.0,observed = voltage_data[0])] 
    for i in range(1,T): 
     V_state.append(Deterministic('V' + str(i), HH(i-1,V_state[i-1],n_state[i-1],m_state[i-1],h_state[i-1])[0])) 
     n_state.append(Deterministic('V' + str(i), HH(i-1,V_state[i-1],n_state[i-1],m_state[i-1],h_state[i-1])[1])) 
     m_state.append(Deterministic('V' + str(i), HH(i-1,V_state[i-1],n_state[i-1],m_state[i-1],h_state[i-1])[2])) 
     h_state.append(Deterministic('V' + str(i), HH(i-1,V_state[i-1],n_state[i-1],m_state[i-1],h_state[i-1])[3])) 
     #A.append(Normal('A' + str(i),mu=V_state[i],sd=2.0,observed = voltage_data[0])) 
    A = Normal('A',mu=V_state,sd=2.0,observed = voltage_data) 
    start = find_MAP() 
    #step = NUTS(scaling=start) 
    step = Metropolis(start=start) 
    trace = sample(100,step) 

是我見過問這些論壇上唯一的另一個例子是在這裏：Simple Dynamical Model in PyMC3

然而，這並不能幫助回答我的問題，因爲無論提出的解決方案有作品。我的解決方案也行不通 - 我沒有收到錯誤，但是當我運行腳本時，抽樣似乎永遠不會開始。無論如何，我的解決方案似乎效率低下，因爲我必須運行python for loop來定義所有的Deterministic分佈。我不知道，如果pymc3甚至認識到我把它們全部放在一個純Python列表中的方式。如果我的函數HH（）可以返回一個numpy數組或某種類型的theano對象，也許這會有所幫助。但是我失去了如何實現它。

Answer 1

你的方法也適用於我，在最後兩行一個微小的變化經過閱讀：

step = Metropolis() 
trace = sample(100, step, start=start) 

這需要一段時間，不過，這樣就可以使T較小，如果你想看到的結果越快。

這裏是a notebook that shows it in action。