functions{ //ODE system vector logistic(real t, //time vector y, // state real r, //average growth rate real K, //carrying capacity / maximum total density of both y real s1, // relative growth (dis)advantage of y1 real s2 // relative growth (dis)advantage of y2 ) { vector[2] dydt; dydt[1] = r*(1+s1)*y[1]*(1-(y[1]+y[2])/K); dydt[2] = r*(1+s2)*y[2]*(1-(y[1]+y[2])/K); return dydt; } // function for splitting ODE system by experiment and fitting real partial_sum(array[] vector y_slice, int start, int end, int nrowz, vector y0,real t0,array[] real ts,real r,real K, array[] real s1, array[] real s2, array[] int IDX, real sigma){ real which_s1 = s1[IDX[end]]; // subsetting using 'end' value of reduce_sum though all values start:end should be identical for any IDX real which_s2 = s2[IDX[end]]; array[nrowz] vector[2] mu = ode_rk45(logistic,y0,t0,ts[start:end],r,K,which_s1,which_s2); //numerical integration step real likely = 0; // initializing likelihood likely += normal_lpdf(y_slice[,1]|mu[,1],sigma); //likelihood for competitor y1 likely += normal_lpdf(y_slice[,2]|mu[,2],sigma); //likelikhood for competitor y2 return likely; } } data{ int N; //rows in data int N_ID; //number of experiments array[N] int IDX; //experiment ID real t0; //initial timepoint (0) array[N] vector[2] yobs; //ob array[N] real ts; // timepoints } parameters{ real inoculum; //y0 is the same for both competitors (experimental constraint) real r; // average growth rate of all competition experiments real K; // carrying capacity / max density array[N_ID] real s1; //array of growth (dis)advantages for all the different y1 tested array[N_ID] real s2; real sigma; //experimental noise hyperparameter } model{ //priors sigma ~ exponential(1); r ~ normal(1,1); //mostly 0.5-1.5 K ~ normal(1,1); //mostly 0.5-1.5 s1 ~ normal(0,0.1); // expecting small differences from mean growth rate s2 ~ normal(0,0.1); inoculum ~ normal(0.03,0.03); //initialize y0 parameter constrained to be same for both competitors vector[2] y0 = rep_vector(inoculum,2); //manual calculation of grain size (all experiments have the same number of timepoints thankfully) int grainsize = N%/%N_ID; //likelihood reduce_sum loop target += reduce_sum_static(partial_sum,yobs,grainsize, grainsize,//second time specifying grainsize is to get number of rows //(is that my problem: I need grainsize to be fixed rather than upper limit?) y0,t0,ts,r,K,s1,s2,IDX,sigma); }