Hello! I am trying to simulate a 3-state Markov chain and the transition probabilities come from a Gaussian process with spatio-temporal covariance matrix and some other variables. I face these two different errors and therefore it says that sampling can not be done.
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: Exception: matrix[uni] indexing: accessing element out of range. index 1 out of range; expecting index to be between 1 and -199758416 (in 'string', line 125, column 8 to column 77)
[1] "Error : Exception: matrix[uni] indexing: accessing element out of range. index 1 out of range; expecting index to be between 1 and -199758416 (in 'string', line 125, column 8 to column 77)"
[1] "error occurred during calling the sampler; sampling not done"
and sometimes this error:
Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: Exception: std::bad_alloc (in 'string', line 129, column 8 to column 77) [origin: bad_alloc]
[1] "Error : Exception: std::bad_alloc (in 'string', line 129, column 8 to column 77) [origin: bad_alloc]"
[1] "error occurred during calling the sampler; sampling not done"
Maybe the 2nd error i because the element is out of range? I am not sure! Any help is appreciated specially on how I can change my code so that I don’t have any memory issues.
Here is the code:
functions {
array[] matrix cor_sep(real nu, real c, real a, real alpha,
data array[] matrix u_ar, data array[] matrix h_ar, data int n_var, data int n_block_row) {
array[n_block_row] matrix[n_var, n_var] corr;
for (n in 1:n_block_row) {
corr[n] = (a * u_ar[n] ^ (2 * alpha) + 1) ^ (-1) .* ((1 - nu) * exp(-c * h_ar[n]));
for (j in 1:n_var) {
for (i in 1:n_var) {
if (i == j)
corr[n][i, j] = (a * (u_ar[n][i, j]) ^ (2 * alpha) + 1) ^ (-1);
}
}
}
return corr;
}
matrix cor_joint(real nu, real c, real a, real alpha, array[] matrix corr_ar,
data int n_var, data int n_block_row) {
int ind_i_start;
int ind_i_end;
int ind_j_start;
int ind_j_end;
matrix[n_var * n_block_row, n_var * n_block_row] corr;
for (i in 1:n_block_row) {
ind_i_start = ((i - 1) * n_var + 1);
ind_i_end = (n_var * i);
for (j in 1:n_block_row) {
ind_j_start = ((j - 1) * n_var + 1);
ind_j_end = (n_var * j);
if (j >= i)
corr[ind_i_start:ind_i_end, ind_j_start:ind_j_end] = corr_ar[j - i + 1];
else
corr[ind_i_start:ind_i_end, ind_j_start:ind_j_end] = corr_ar[i - j + 1];
}
}
return corr;
}
}
data {
int<lower=1> N;
int<lower=1> N_mu;
int<lower=1> lag_u;
int<lower=1> n_ahead;
//number of wind farms
int<lower=1> n_var;
int<lower=1> lag_max;
int<lower=1> n_block_row;
matrix[n_var, n_var] dist;
vector[N_mu] X[N];
array[n_block_row] matrix[n_var, n_var] h_ar;
array[n_block_row] matrix[n_var, n_var] u_ar;
//the final y
int<lower = 1, upper = 3> y[N,n_var];
}
parameters {
// model parameters
real<lower=0, upper=1> nu;
real<lower=0> c;
real<lower=0> a;
real<lower=0, upper=1> alpha;
//new_parameters
vector[N]beta[3,3];
array[3,3] matrix[N,n_var] b;
// mean
vector[N_mu] mu;
}
transformed parameters {
array[n_block_row] matrix[n_var, n_var] c_sep;
matrix[n_var * n_block_row, n_var * n_block_row] c_joint;
array[3,3] matrix[N,n_var] theta;
c_sep = cor_sep(nu, c, a, alpha, u_ar, h_ar, n_var, n_block_row);
c_joint = cor_joint(nu, c, a, alpha, c_sep, n_var, n_block_row);
for(t in 1:3){
for(k in 1:3){
for(i in 1:N){
for(j in 1:n_var){
theta[t,k,i,j] = beta[t,k,i]*X[i,j] + b[t,k,i,j];
}
}
}
}
}
model {
// flat priors
nu ~ uniform(0, 1);
c ~ uniform(0.0001, 5);
a ~ uniform(0.0001, 5);
alpha ~ uniform(0, 1);
mu ~ uniform(-10,10);
for(t in 1:3){
for(k in 1:3){
for(i in 2:N){
for(j in 1:n_var){
b[t,k,i,j]~ uniform(0.001, 1);
beta[t,k,i] ~ uniform(0.001, 1);
}
}
}
}
for(i in 2:N){
X[i]~ multi_normal(mu, c_joint);
for(j in 1:n_var){
y[i,j] ~ categorical(softmax(to_vector(theta[y[i-1,j]][1:3][i][j])));;
}
}
}