Hi, I am trying to run MCMC sampling with Turing.jl for a system of ODEs in which some parameters depend on an external variable (voltage). Furthermore, I would like to use data from multiple different types of experiments, which correspond to different outputs of the model.
However, it seems that the functions I use to get the outputs I need from my model are not compatible with Turing.jl, but the error messages overflow my terminal, so I haven’t been able to identify specific problems. Any advice would be greatly appreciated!
Briefly, I have four state variables (A, B, C, D), and my data reflect the time evolution of the sum of C and D. The different types of data include: (1) responses of the system to different individual voltages when initialized from one fixed voltage, and (2) responses of the system to more complex patterns of voltage.
Focusing only on the former type of data, my approach right now (which doesn’t work for MCMC) was to add voltage as a parameter, solve an ODEProblem for each voltage, and collect the sum C + D in an array:
#a system of ODEs
function f!(du, u, p, t)
V, a, b, c, ... = p
a1 = a*exp(-V)
b1 = b*exp(V)
...
mul!(du,A,u) #A is a matrix containing appropriate coefficients for the ODEs)
end
p_true = [...] #parameters
trange = 1000 #simulation time
vrange = 0:10:100 #simulation voltages
tspan = (1.0, length(trange)) #simulation time
#copy parameters and insert 0 in the first position to act as a placeholder for voltage later on
p_true2 = copy(p_true)
pushfirst!(p_true2, 0)
# dt = number of timepoints, dv = number of voltages
my_data = zeros((dt, dv))
for i = 1
j = 1
x0 = get_ss(p_true, initial_voltage) #steady state at the initial voltage
while j <= dv
p_true2[1] = vrange[j] #update voltage in-place
prob = ODEProblem(f!, x0, tspan, p_true2, saveat=1.0)
sol = solve(prob)
my_data[:,j] = sol[3,:] + sol[4,:] #get the sum C+D
j += 1
end
end
This works pretty well to just run the model. For MCMC, I basically copied the above into a Turing.jl model, e.g.
@model function fit(data)
σ ~ InverseGamma(2, 3)
a ~ pd[1] #pd is a list of Distributions
b ~ pd[2]
...
p = [a, b, c, ...] #sampled parameters
p2 = copy(p) #do the same thing as above, i.e. create a placeholder for voltage in the first position of the parameters vector
pushfirst!(p2, 0)
predicted = (collection of C+D sums after solving an ODEProblem for each voltage, as above)
for j = 1:size(predicted)[2]
for i = 1:size(predicted)[1]
data[i,j] ~ Normal(predicted[i,j], σ)
end
end
end
When I try to run this, I think I’m getting an error related to autodifferentiation, but I’m not sure - maybe I should be using a different tool/package?
I have been able to reproduce some simpler examples such as the Lotka-Volterra equations, but I couldn’t find any examples that were helpful for my use case above, where I believe each voltage yields a different set of parameters and therefore a different dataset (?).