Réécriture d'un script pymc pour l'estimation de paramètres dans les systèmes dynamiques de pymc3

Question

Je voudrais utiliser pymc3 pour estimer des paramètres et des états inconnus dans un modèle de neurones de Hodgkin Huxley. Mon code dans pymc est basé sur http://healthyalgorithms.com/2010/10/19/mcmc-in-python-how-to-stick-a-statistical-model-on-a-system-dynamics-model-in-pymc/ et s'exécute raisonnablement bien.

#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... 

Dans pymc3, l'astuce Lambda n'est pas disponible. Voici une de mes tentatives:

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) 

Le seul autre exemple que j'ai vu demandé sur ces forums est ici: Simple Dynamical Model in PyMC3

Cependant, cela ne permet pas de répondre à ma question parce qu'il ya la solution proposée ne travaux. Ma solution ne fonctionne pas non plus - je n'ai pas d'erreur, mais quand je lance le script, l'échantillonnage ne semble jamais commencer. Dans tous les cas, ma solution semble inefficace car je dois lancer une boucle python pour définir toutes les distributions déterministes. Je ne suis pas sûr si pymc3 reconnaît même la façon dont je les ai tous mis dans une liste python pure. Si ma fonction HH() pouvait renvoyer un tableau numpy ou une sorte d'objet theano, peut-être que cela aiderait. Mais je suis perdu quant à la façon de l'implémenter.

Answer 1

Votre approche fonctionne pour moi, après un petit changement dans les deux dernières lignes à lire:

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

Il prend un certain temps, cependant, pour que vous puissiez faire T plus petit si vous voulez voir des résultats plus tôt.

Voici a notebook that shows it in action.