Решение ODE в PYMC3
Здесь я хочу оценить параметры (гамма и омега) затухающего гармонического осциллятора, заданного
йХ ^2/ дт ^2+ гамма * аЯ / дт +(2* пи * омега)^2* Х =0. (Мы можем добавить белый гауссовский шум в систему.)
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_model.py, мы можем запустить PYMC следующим образом
import pymc
import DampedOscil_model
MDL = pymc.MCMC(DampedOscil_model, db='pickle')
MDL.sample(iter=1e4, burn=1e2, thin=2)
gama_trace=MDL.trace('gama')[- 1000:]
omega_trace=MDL.trace('omega')[-1000:]
gama=MDL.gama.value
omega=MDL.omega.value
И это хорошо работает (см. Ниже).
Истинный сигнал, построенный по gama_true=2.0 и omega_est=1.5 по сравнению с оцененным сигналом. Расчетные значения параметров: gama_est=2,04 и 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.Deterministic('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 нет лямбда-функции.
Я ценю вашу помощь в решении ODEs в PYMC3 для решения задачи оценки параметров в биологических системах (оценка параметров уравнения по данным).
Заранее большое спасибо за вашу помощь.
С уважением,
Meysam
0 ответов
Я бы предложил и успешно реализовал использование метода "черного ящика" для взаимодействия с PYMC3. В данном случае это означает самостоятельное вычисление логарифмической вероятности с последующим использованием PYMC3 для ее выборки. Это требует написания ваших функций таким образом, чтобы Theano и PYMC3 могли взаимодействовать с ними.
Это описано в блокноте на странице PYMC3, где в качестве примера используется cython.
Вот немного более короткий пример того, что нужно сделать.
Сначала вы можете загрузить свои данные и настроить любые необходимые параметры, такие как ваши временные шаги и т. Д.
import pymc3 as pm
import numpy as np
import theano
import theano.tensor as tt
#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]
Затем вы определяете функцию генерации данных так же, как и раньше, но вам не нужно использовать для этого какие-либо декораторы из PYMC. Результатом этой функции должно быть то, что вам нужно сравнить с вашими данными, чтобы рассчитать вероятность.
def DHOS(theta):
gama,omega=theta
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]
Затем вы пишете функцию, которая вызывает предыдущую функцию и вычисляет вероятность, используя любое распределение, которое вы хотите, в этом нормальном распределении.
def my_loglike(theta,data,sigma):
"""
A Gaussian log-likelihood function for a model with parameters given in theta
"""
model = DHOS(theta) #V1 and V2 from the DHOS function
#Here data = [ydata1,ydata2] to compare with model
#sigma is either the same shape as model or a scalar
#which corresponds to the uncertainty on the data.
return -(0.5)*sum((data - model)**2/sigma**2)
Отсюда вам нужно определить класс Theano, чтобы он мог взаимодействовать с PYMC3.
# define a theano Op for our likelihood function
class LogLike(tt.Op):
"""
Specify what type of object will be passed and returned to the Op when it is
called. In our case we will be passing it a vector of values (the parameters
that define our model) and returning a single "scalar" value (the
log-likelihood)
"""
itypes = [tt.dvector] # expects a vector of parameter values when called
otypes = [tt.dscalar] # outputs a single scalar value (the log likelihood)
def __init__(self, loglike, data, sigma):
"""
Initialise the Op with various things that our log-likelihood function
requires. Below are the things that are needed in this particular
example.
Parameters
----------
loglike:
The log-likelihood (or whatever) function we've defined
data:
The "observed" data that our log-likelihood function takes in
x:
The dependent variable (aka 'x') that our model requires
sigma:
The noise standard deviation that our function requires.
"""
# add inputs as class attributes
self.likelihood = loglike
self.data = data
self.sigma = sigma
def perform(self, node, inputs, outputs):
# the method that is used when calling the Op
theta, = inputs # this will contain my variables
# call the log-likelihood function
logl = self.likelihood(theta, self.data, self.sigma)
outputs[0][0] = array(logl) # output the log-likelihood
Наконец, вы можете использовать PYMC3 для построения вашей модели и соответствующего образца.
ndraws = 10000 # number of draws from the distribution
nburn = 1000 # number of "burn-in points" (which we'll discard)
# create our Op
logl = LogLike(my_loglike, rdat_sim, 10)
# use PyMC3 to sampler from log-likelihood
with pm.Model():
gama= pm.Uniform('gama', 0.0, 20.0)
omega=pm.Uniform('omega',0.0, 20.0)
# convert m and c to a tensor vector
theta = tt.as_tensor_variable([gama,omega])
# use a DensityDist (use a lamdba function to "call" the Op)
pm.DensityDist('likelihood', lambda v: logl(v), observed={'v': theta})
trace = pm.sample(ndraws, tune=nburn, discard_tuned_samples=True)
И вы можете использовать внутреннее построение, чтобы увидеть результаты выборки
_ = pm.traceplot(trace)
This was just adapted from the example notebook in the link, and as mentioned there if you want to use NUTS you need gradient information, which you do not have given you custom function. In the link it talks about how to sample the gradient and construct it so you can pass it into the sampler, but I have not shown that here.
Additionally if you want to use solve_ivp (or odeint or another solver), all you have to do is change the DHOS function as you normally would to invoke the solver. The rest of the code should be portable to whatever problem you, or anyone else, need.