functions{ matrix vcovFBM(vector t,//input timeSince_t0 from ni_l to ni_u real H, int n,//input n as an iteger real kappa) { matrix[n, n] m; for(j in 1:n) { for(k in 1:n) { if(j == k) { m[j, k] = pow(fabs(t[j]), 2*H); } else{ m[j, k] = kappa*0.5*(pow(fabs(t[j]), 2*H) + pow(fabs(t[k]), 2*H) - pow(fabs(t[j] - t[k]), 2*H)); } } } return (cholesky_decompose(m)); } } data { // Number of level-1 observations (an integer) int N; // Number of level-2 clusters int Ni; // Cluster IDs int level2[N]; //Number of predictors int K; // Continuous outcome real Y_ij[N]; // Continuous predictor matrix[N, K] X_ij; //used for covFBM int ni_l[Ni]; int ni_u[Ni]; vector[N] timeSince_t0; int ni[Ni]; } parameters { // Define parameters to estimate //Parameters of covariates vector[K] alpha; // Level-1 real sigma; // Level-2 random effect real ui[Ni]; real tau; vector[N] wij; real kappa; real H; } transformed parameters { // Individual mean vector[N] mu; // Individual mean for (i in 1:N) { mu[i] = ui[level2[i]] + X_ij[i,]*alpha + wij[i]; } } model { // Prior part of Bayesian inference sigma ~ inv_gamma(1, 10); tau ~ inv_gamma(1, 10); alpha ~ normal(0, 10); H ~ uniform(0, 1); kappa ~ uniform(0, 1); // Random effects distribution ui ~ normal(0, tau); //FBM for(i in 1 : Ni) { wij[ni_l[i]:ni_u[i]] ~ multi_normal_cholesky(rep_vector(0, ni[i]), vcovFBM(timeSince_t0[ni_l[i]:ni_u[i]], H, ni[i], kappa)); } // Likelihood part of Bayesian inference // Outcome model N(mu, sigma^2) (use SD rather than Var) for (i in 1:N) { Y_ij[i] ~ normal_lpdf(mu[i], sigma); } }