Help with a dynamic system

Hi all, I’m trying to get familiar with Stan, but have difficulties building a model. I looked in different topics and documentations but couldn’t find a template, which fits my needs.

What I want to do: I have a dynamic system, consisting of ca. 200 differential equations, containing 50 parameters I need to sample. The measure for the data is a constellation of some states at different timepoints for 10 different initial states. I hope I can explain it with an example: I know in the experiment with initial state [1,2,1,1,3,4,…,1] that at t=50s I measured x[4]+x[79]+x[97]+x[179]=153. There is no way to measure the states separately. I know from the nature of the experiment that my measure is normal distributed, where the sigma and the scale is unknown. So it would be great to have also a sampling for sigma and scale.

I tried to build the model, but can’t figure out, how to implement the data. Here is my attempt with a smaller artificial model.

functions {
	real[] 	reactions(real t,
				real[] x,
				real[] theta) {
		real dxdt[10];
		
		dxdt[1] =  theta[1]*x[2] - theta[2]*x[1];
		dxdt[2] = - theta[3]*x[2];
		dxdt[3] =  theta[4]*x[4]*x[4] - theta[5]*x[1];
		dxdt[4] = - theta[2]*x[4];
                dxdt[5] =  theta[4]*x[5] - theta[6]*x[2];
		dxdt[6] = - theta[1]*x[8];
                dxdt[7] =  theta[3]*x[6] - theta[6]*x[2];
		dxdt[8] =  theta[5]*x[10];
                dxdt[9] =  theta[3]*x[8] + theta[2]*x[9];
		dxdt[10] =  theta[6]*x[9];
		return dxdt;
	}
}
data {
	int<lower=1> T;
	int<lower=1> C;
	real<lower=0> measure[T,C,3];
        real initial[C,10]
	real ts[T];
}

parameters {
	vector<lower=0>[8] sigma;
	real<lower=-5,upper=5> theta[6];
}
transformed parameters {
	real x_hat[T,2];	
        for(c in 1:C)
             x_hat[c] = integrate_ode_bdf(reactions, initial[c], 0, ts, theta);
}
model{
	sigma~cauchy(0,1);
	for (t in 1:T)
            for(c in 1:C)
                measure[t,c,1]~normal(x_hat[c,2]+x_hat[c,5]+x_hat[c,6],sigma);
                measure[t,c,2]~normal(x_hat[c,1]+x_hat[c,4]+x_hat[c,6],sigma);
                measure[t,c,3]~normal(x_hat[c,4]+x_hat[c,7],sigma);
}

Obviously this doesn’t work at all, but I would be happy for any help you can provide. This doesn’t help me: http://mc-stan.org/events/stancon2017-notebooks/stancon2017-margossian-gillespie-ode.html .

That seems like the right general approach. We don’t have efficient ways to avoid Jacobian calculations for the x_hat you don’t need, but that shouldn’t stop this from working.

Could you say more about what kinds of problems you’re having? Have you managed to get any kind of noisy diff eq model to work in Stan? It might help to build up from easier pieces to see how everything fits together.