Hi,
I am trying to build up the following multivariate hierarchical model using an ODE. The data set consists of time series measurements for J subjects. Thus for each subject observed measurements (y_{j}) are modeled as:
y_j = f_{j}(\theta_{j}) + e_{j}, \quad j = 1:J, \theta_{j} = (\theta_{1}, \dots, \theta_{K}) ,
where f_{j}(\cdot ) is solution of ODE for parameter set \theta_{j} and e_{j} \sim MultiNormal(0, \sigma^{2} I{n_{j}}), with I as identity matrix and n_{j} is the number of measurements for subject j.
And log(\theta_{j}) = a +b_{j}, where b_{j} \sim MultiNormal(0, \Sigma), a \sim MultiNormal(mu_a, \Lambda) with \Lambda a diagonal matrix. (Note: I have logged parameters as some parameter values are large and some are very small).
At this stage, I am coding the priors and computing log(\theta_{j}) . Following Stans manual, I gave \Sigma an LKJ prior. The code runs okay, but I am stuck with constraining one parameter (parameter d \in [0,1]) in the parameter set, which must be between 0 and 1 on the original (non logged) scale.
Currently, taking exponential of l_theta give values for d larger than 1. How should I modify the code so that d is between 0 and 1 on the original scale?
functions {
real[] ode(real t,
real[] y,
real[] theta,
real[] x_r, // x_r:real data, x_i: integer data
int[] x_i) {
real dydt[4];
real A = theta[1];
real beta = theta[2];
real delta = theta[3];
real gamma = theta [4];
real omega = theta[5];
real alpha = theta[6];
real kappa = theta[7];
real eta = theta[8];
real d = theta[9];
dydt[1] = A - gamma*y[1]-beta * y[3] * y[1];
dydt[2] = beta * y[1] * y[3] - delta * y[2];
dydt[3] = (1-d)*omega * y[2] - kappa * y[3] -alpha* y[4]*y[3] ;
dydt[4] = eta *y[3]* y[4];
return dydt;
}
}
data {
int<lower = 1> J; //number of subjects
int<lower =1> K; //number of parameters
vector[K] mu_a; //mean vector of a
}
transformed data {
real x_r[0];
int x_i[0];
}
parameters{
matrix[J,K] a;
vector<lower = 0>[K] tau; // sd of b
cholesky_factor_corr [K] Omega;
matrix[K,J] z;
}
transformed parameters {
matrix[J,K] b;
matrix[J,K] l_theta; //parameters on the log scale
b = (diag_pre_multiply(tau, Omega)*z)';
l_theta = a + b;
}
model {
//priors
tau ~ cauchy(0,2.5);
Omega ~ lkj_corr_cholesky(2);
to_vector(z) ~ normal(0,1);//this implies for(k in 1:K) z[k]~normal(0,1)
to_vector(a) ~ normal(mu_a, 5);
}
The R code used to run the stan program is :
stan_data <- list( J = 1, #numb. of subjects
K = 9, #numb. of parameters
mu_a = c(log(1.4e7), log(1.72e-10), log(0.14),log(0.14), log(1e4), log(0.001),log(3.48),log(1.29e-5), log(0.5))
#mean vector of (A, beta, delta, gamma,omega,alpha, kappa,eta, d) on log scale
)
fake_data <- stan("model_fit.stan",
data = stan_data,
chains = 1, iter = 1000,
seed=4838282)
a <- extract(fake_data, pars =c("a", 'b', 'l_theta'))
Also, is the following code used to sample a \sim MultiNormal(mu_a , \Lambda), with \Lambda = diag(rep(10, K)) correct? I read from Stan manual to use a normal distribution if the co-variance matrix of multinormal is diagonal.
to_vector(a) ~ normal(mu_a, 5);
I feel that there is something incorrect with the code as the density of l_theta for each parameter is concentrated at 0 and taking exponential do not give a proper distribution for theta parameters.
plot(density(exp(a$l_theta[,1,7])))
sampling
Thanks