Hi,
I’m trying to model activity propagation in a network with N nodes for N_TS time points that spreads as described by the following evolution equation:
\frac{dx_i}{dt}\tau_i = -x_i(t) + s\phi\big(x_i(t)\big) + g\Bigg(\sum_{j\neq i}^{N} W_{ij}\phi\big(x_j(t)\big)\Bigg) + I_i(t)
This is the first time I’m trying to solve an ODE in a Stan model so it’s very possible that I’ve made other mistakes as well. I haven’t even gotten to the point of testing the actual model yet because I can’t compile it due to the following error:
Compiling Stan program...
/
Semantic error in '/var/folders/10/zpzbg__x3hlc56gwq9zf6r600000gn/T/RtmphjhKZn/model-e51f176820ee.stan', line 120, column 22 to column 82:
-------------------------------------------------
118: // states x are estimated as parameters z with some uncertainty.
119: // these estimated parameters are used in the model description to relate them to measured data
120: vector[N] z[N_TS] = ode_rk45(dx_dt, z_init, t0, ts, N, N_t, I, ts, s, g, b, tau);
^
121: }
122:
-------------------------------------------------
Ill-typed arguments supplied to function 'ode_rk45'. Available signatures:
(real, vector, ...) => vector, vector, real, real[], ...
.
Instead supplied arguments of incompatible type:
(real, vector, int, vector[], vector, vector, real, real, real, real) => vector, vector, real, real[], int, vector[], vector, real[], real, real, real, real.
The supplied arguments of type real and vector look like they should agree with the required signature for the ode_rk45 function so I’m not sure how to debug this. As I said, it is entirely possible that the error lies somewhere else but I’m not sure where to start. Any help would be appreciated.
Here’s the full model:
functions {
real phi(real a){
real out = (exp(2*a)-1)/(exp(2*a)+1);
return out;
}
// From
// https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/1X/25b0af2e6aba491a267c45fa090c6374ad57c4b4.pdf
// https://discourse.mc-stan.org/t/zero-chain-movement-mixing-in-forced-ode-model/921
int find_interval_elem(real x, vector sorted, int start_ind){
int res;
int N;
int max_iter;
real left;
real right;
int left_ind;
int right_ind;
int iter;
N = num_elements(sorted);
if(N == 0) return(0);
left_ind = start_ind;
right_ind = N;
max_iter = 100 * N;
left = sorted[left_ind ] - x;
right = sorted[right_ind] - x;
if(0 <= left) return(left_ind-1);
if(0 == right) return(N-1);
if(0 > right) return(N);
iter = 1;
while((right_ind - left_ind) > 1 && iter != max_iter) {
int mid_ind;
real mid;
// is there a controlled way without being yelled at with a
// warning?
mid_ind = (left_ind + right_ind) / 2;
mid = sorted[mid_ind] - x;
if (mid == 0) return(mid_ind-1);
if (left * mid < 0) { right = mid; right_ind = mid_ind; }
if (right * mid < 0) { left = mid; left_ind = mid_ind; }
iter = iter + 1;
}
if(iter == max_iter)
print("Maximum number of iterations reached.");
return(left_ind);
}
vector dx_dt(real t, // time
vector x, // system state {1 x N - network activity for each time point}
int N,
vector[] N_t, //vector array
vector I,
vector ts,
real s,
real g,
real b,
real tau) {
vector[N] res;
int d = find_interval_elem(t, ts, 1);
for (i in 1:N){
res[i] = (-x[i] + s * phi(x[i]) + g *N_t[d, i] + b * I[d])/tau;
}
// returns change for one time point for all nodes
return res;
}
}
data {
int<lower = 0> N_TS; // number of measurement times
int<lower = 0> N; // number of nodes
real t0;
real ts[N_TS]; // measurement times > 0
vector[N] y_init; // initial measured activity level
// Arrays are declared by enclosing the dimensions in square brackets following the name of the variable.
// this is an array of length N_TS consisting of (column) vectors of length N
// Each array element is like a row
vector[N] y[N_TS]; // measured activity level with nodes in cols and timepoints in rows
matrix[N, N] W;// adjacency matrix
vector[N_TS] I; // task stimulation
}
transformed data{
vector[N] N_t[N_TS];
for(i in 1:N){
for(t in 1:N_TS){
//row of incoming connections to node i of length n
//multiplied with activity (column) vector of length n for all nodes at time t
N_t[t, i] = dot_product(W[i,], y[t]);
}
}
}
parameters {
real<lower = 0> s; // self coupling
real<lower = 0> g; // global coupling
real<lower = 0> b; // task modulation
real<lower = 0> tau; // time constant
real<lower = 0> sigma; // measurement error
vector[N] x_init;
}
transformed parameters {
// states x are estimated as parameters z with some uncertainty.
// these estimated parameters are used in the model description to relate them to measured data
vector[N] x[N_TS] = ode_rk45(dx_dt, x_init, t0, ts, N, N_t, I, ts, s, g, b, tau);
}
model {
s ~ lognormal(-1, 1);
g ~ lognormal(-1, 1);
beta ~ lognormal(-1, 1);
tau ~ normal(1, 0.5);
sigma ~ lognormal(-1, 1);
for (k in 1:N) {
y_init[k] ~ lognormal(log(x_init[k]), sigma);
y[ , k] ~ lognormal(log(x[, k]), sigma[k]);
}
}