python – 在pymc3中重写动态系统中参数估计的pymc脚本

概述我想使用pymc3来估计霍奇金赫胥黎神经元模型中的未知参数和状态.我在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@de 我想使用pymc3来估计霍奇金赫胥黎神经元模型中的未知参数和状态.我在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@deterministicdef 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 likelihoodA = normal('A',mu=V,tau=2.0,value=voltage_data,observed = True)#run the sampling...

在pymc3中,Lambda技巧不适用于我.这是我的一个尝试：

import numpy as npimport theano.tensor as ttfrom pymc3 import Model,normal,Uniform,Deterministic,sample,NUTS,Metropolis,find_MAPimport scipy#observed dataT = 10dt = 0.02voltage_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 laterdef 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))#ModelNaKL_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',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',-30.0)    dvm = Uniform('dvm',10.0,20.0)    dvm_t = Uniform('dvm_t',20.0)    tmm_0 = Uniform('tmm_0',0.05,0.25)    tm_1 = Uniform('tm_1',1.0)    vh = Uniform('vh',-80.0,-40.0)    vh_t = Uniform('vh_t',-40.0)    dvh = Uniform('dvh',-30.0,-10.0)    dvh_t = Uniform('dvh_t',-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',-40.0)    dvn = Uniform('dvn',40.0)    dvn_t = Uniform('dvn_t',40.0)    tn_0 = Uniform('tn_0',1.5)    tn_1 = Uniform('tn_1',3.0,7.0)    #Initial Conditions    V_0 = Uniform('V_0',50.0)    n_0 = Uniform('n_0',1.0)    m_0 = Uniform('m_0',1.0)    h_0 = Uniform('h_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,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),h_state[i-1])[1]))        m_state.append(Deterministic('V' + str(i),h_state[i-1])[2]))        h_state.append(Deterministic('V' + str(i),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,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循环来定义所有的Deterministic发行版.我不确定pymc3是否能够识别我将它们全部放在纯python列表中的方式.如果我的函数HH()可以返回一个numpy数组或某种theano对象,也许这会有所帮助.但我迷失了如何实现它.

解决方法 在最后两行的微小变化后,您的方法对我有用：

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

但是,这需要一段时间,因此如果您希望更快地看到结果,可以使T更小.

这是a notebook that shows it in action.

总结

以上是内存溢出为你收集整理的python – 在pymc3中重写动态系统中参数估计的pymc脚本全部内容，希望文章能够帮你解决python – 在pymc3中重写动态系统中参数估计的pymc脚本所遇到的程序开发问题。

如果觉得内存溢出网站内容还不错，欢迎将内存溢出网站推荐给程序员好友。