在PYMC3中求解ODEs

Question

这里我的目的是估计阻尼谐振子的参数（gama和omega）由

dX^2/dt^2+gamma*dX/dt+(2*pi*omega)^2*X=0. (We can add white gaussian noise to the system.)

---

 import pymc
 import numpy as np
 import scipy.io as sio
 import matplotlib.pyplot as plt; 
 from scipy.integrate import odeint

 #import data
 xdata = sio.loadmat('T.mat')['T'][0]  #time
 ydata1 = sio.loadmat('V1.mat')['V1'][0]  #  V2=dV1/dt, (X=V1),
 ydata2 = sio.loadmat('V2.mat')['V2'][0]  # dV2/dt=-(2pi*omega)^2*V1-gama*V2

 #time span for solving the equations
 npts= 500   
 dt=0.01
 Tspan=5.0
 time = np.linspace(0,Tspan,npts+1) 

 #initial condition
 V0 = [1.0, 1.0]

# Priors for unknown model parameters
sigma = pymc.Uniform('sigma',  0.0, 100.0)
gama= pymc.Uniform('gama', 0.0, 20.0) 
omega=pymc.Uniform('omega',0.0, 20.0)

#Solve the equations
@pymc.deterministic
def DHOS(gama=gama, omega=omega):
V1= np.zeros(npts+1)
V2= np.zeros(npts+1)
V1[0] = V0[0]
V2[0] = V0[1]
for i in range(1,npts+1):
    V1[i]=  V1[i-1] + dt*V2[i-1]; 
    V2[i] = V2[i-1] + dt*(-((2*np.pi*omega)**2)*V1[i-1]-gama*V2[i-1]); 
return [V1, V2]

#or we can use odeint
#@pymc.deterministic
#def DHS( gama=gama, omega=omega):
#    def DOS_func(y, time):
#        V1, V2 = y[0], y[1]
#        dV1dt = V2
#        dV2dt = -((2*np.pi*omega)**2)* V1 -gama*V2 
#        dydt = [dV1dt, dV2dt]
#        return dydt
#    soln = odeint(DOS_func,V0, time)
#    V1, V2 = soln[:,0], soln[:,1]
#    return V1, V2


#  value of outcome (observations)
V1 = pymc.Lambda('V1', lambda DHOS=DHOS: DHOS[0])
V2 = pymc.Lambda('V2', lambda DHOS=DHOS: DHOS[1])


# liklihood of observations
Yobs1 = pymc.Normal('Yobs1', mu=V1, tau=1.0/sigma**2, value=ydata1, observed=True)
Yobs2 = pymc.Normal('Yobs2', mu=V2, tau=1.0/sigma**2, value=ydata2, observed=True)

通过将上述代码保存为DampedOscil_模型.py，然后我们可以运行PYMC，如下所示

^{pr2}$

---

而且效果很好（见下文）。在

The true signal constructed by gama_true=2.0 and omega_est=1.5 versus the estimated signal. The estimated parameter values are gama_est=2.04 and omega_est=1.49

现在我将把这个代码转换成PYMC3来使用NUTS和ADVI。在

import matplotlib.pyplot as plt
import scipy.io as sio
import pandas as pd
import numpy as np

import pymc3 as pm
import theano.tensor as tt
import theano

from pymc3 import Model, Normal, HalfNormal, Uniform
from pymc3 import NUTS, find_MAP, sample, Slice, traceplot, summary
from pymc3 import Deterministic  
from scipy.optimize import fmin_powell


#import data
xdata = sio.loadmat('T.mat')['T'][0]  #time
ydata1 = sio.loadmat('V1.mat')['V1'][0]  #  V2=dV1/dt, (X=V1),
ydata2 = sio.loadmat('V2.mat')['V2'][0]  # dV2/dt=-(2pi*omega)^2*V1-gama*V2

#time span for solving the equations
npts= 500   
dt=0.01
Tspan=5.0
time = np.linspace(0,Tspan,npts+1) 

niter=10000
burn=niter//2;

with pm.Model() as model:

 #Priors for unknown model parameters
 sigma = pm.HalfNormal('sigma', sd=1)
 gama= pm.Uniform('gama', 0.0, 20.0)
 omega=pm.Uniform('omega',0.0, 20.0)

 #initial condition
 V0 = [1.0, 1.0]

#Solve the equations     
# do I need to use theano.tensor here?!
 @theano.compile.ops.as_op(itypes=[tt.dscalar, tt.dscalar],otypes=[tt.dvector])  
def DHOS(gama=gama, omega=omega):
V1= np.zeros(npts+1)
V2= np.zeros(npts+1)
V1[0] = V0[0]
V2[0] = V0[1]
for i in range(1,npts+1):
    V1[i]=  V1[i-1] + dt*V2[i-1]; 
    V2[i] = V2[i-1] + dt*(-((2*np.pi*1)**2)*V1[i-1]-gama*V2[i-1]); 
return V1,V2


 V1 = pm.Deterministic('V1', DHOS[0])
 V2 = pm.Deterministic('V2', DHOS[1])


 start = pm.find_MAP(fmin=fmin_powell, disp=True)
 step=pm.NUTS()
 trace=pm.sample(niter, step, start=start, progressbar=False)

 traceplot(trace);

 Summary=pm.df_summary(trace[-1000:])

  gama_trace = trace.get_values('gama', burn)
  omega_trace = trace.get_values('omega', burn)

对于此代码，我得到以下错误： V1=pm.确定性（'V1'，DHOS[0]）

      TypeError: 'FromFunctionOp' object does not support indexing

简而言之，我想知道如何将PYMC代码的以下部分转换为PYMC3。在

@pymc.deterministic
def DOS(gama=gama, omega=omega):
V1= np.zeros(npts+1)
V2= np.zeros(npts+1)
V1[0] = V0[0]
V2[0] = V0[1]
for i in range(1,npts+1):
    V1[i]=  V1[i-1] + dt*V2[i-1]; 
    V2[i] = V2[i-1] + dt*(-((2*np.pi*omega)**2)*V1[i-1]-gama*V2[i-1]); 
return [V1, V2]


V1 = pymc.Lambda('V1', lambda DOS=DOS: DOS[0]) 
V2 = pymc.Lambda('V2', lambda DOS=DOS: DOS[1])

问题是，首先，PYMC3中确定性函数的论证与PYMC不同，其次，PYMC3中没有Lambda函数。在

感谢您在PYMC3中求解ODE，以解决生物系统中的参数估计任务（从数据中估计方程参数）。在

提前非常感谢你的帮助。在

谨致问候

蛹