Hi everyone, I run a lengthy stan code and got stuck at getting index for tensors. The function “mapping” takes an integer index and transforms it into a 3D vector T, which is used as index later in the transformed parameter block. The issue is that the type of T, T, T is real not integer so it gets me the error “Index must be of type int or int or must be a range. Instead found type real.” in " transformed parameters {
real<lower=0> r_n[K, D0, S, U0]; // tensor topic
vector T;
for (k in 1:K)
for (i in 1:I) {
T = mapping(i);
r_n[k, T, T, T] = r[k, i];
}
}".
May I ask how I can modify the function mapping to fix this problem? Thanks!!!

functions {
vector mapping(int I){
int L;
int remainder;
vector T;
L = I - 1;
remainder = L%16;
T = L/16; //label for S
T = remainder/4; //label for D
T = remainder%4;	//label for U
return T;
}
}
data {
int<lower = 1> S;	//number of single-base substitution
int<lower = 1> J;	//number of samples
int<lower = 1> K;   //number of sigs
int<lower = 1> D0;	//number of downstream bases at #0
int<lower = 1> U0;   //number of upstream bases at #0
int<lower = 1> D1;	//number of downstream bases at #1
int<lower = 1> U1;   //number of upstream bases at #1
int<lower = 1> I;	//number of tri-nucleotide substitutions

int<lower = 0> X[D1, D0, S, U0, U1, J];	// the data mat motifys by samples

real<lower=0> alpha;   // hyperparam for signature rates
real<lower=0> beta;    //hyperparam for up/downstream bases rates

real<lower=0> gamma0;  // hyperparam for shape parameters
real<lower=0> gamma1;  // hyperparam for shape parameters
}

transformed data {
vector<lower = 0>[I] alpha_array = rep_vector(alpha, I); //Dirichlet params for signatures
vector<lower = 0>[D0] beta_array = rep_vector(beta, D0);
}

parameters {
simplex[U1] u[K, U0, S];	//inferred upstream bases
simplex[D1] d[K, D0, S];	//inferred downstream bases
simplex[I] r[K];	//inferred sigs K by I
}

transformed parameters {
real<lower=0> r_n[K, D0, S, U0]; 	// tensor topic
vector T;
for (k in 1:K)
for (i in 1:I) {
T = mapping(i);
r_n[k, T, T, T] = r[k, i];
}
}

model {
real mutation_rate;
vector T;
int d0;
int s;
int u0;

for (k in 1:K) {
nu[k] ~ inv_gamma(gamma0, gamma1);
mu[k] ~ inv_gamma(delta0, delta1);
r[k] ~ dirichlet(alpha_array);

for (u0 in 1:U0)
for (s in 1:S){
u[u0, s, k] ~ dirichlet(beta_array);
d[u0, s, k] ~ dirichlet(beta_array);
}
}

for (j in 1:J)
for (k in 1:K) {
theta[k, j] ~ gamma(nu[k], nu[k]/mu[k]);
}

for (j in 1:J)
for (d1 in 1:D1)
for (u1 in 1:U1)
for (i in 1:I){
mutation_rate = 0;
T = mapping(i);
d0 = T;
s = T;
u0 = T;
for (k in 1:K){
mutation_rate += u[k, s, u0, u1]*r_n[k, d0, s, u0]*d[k, d0, s, d1]*theta[k, j];
}
X[d1, d0, s, u0, u1, j] ~ poisson(mutation_rate);

}

}



model{
vector[N] mu = alpha+beta*x;
y~normal(mu,sigma);
}

To include mathematical notation in your post put LaTeX syntax between two \$ symbols, e.g.,
p(\theta | y) \propto p(\theta) p(y | \theta).

Don’t forget to add relevant tags to your topic (top right of this form) for application area and/or class of models you work with.

Rather than vector you can use the type array[] int for the return type of mapping(). Similarly, where you say vector T you can use array int T

The advantage of using arrays is that if you have

array[M] vector[N] x;


then the operation x[m] can operate by reference (i.e., it’s constant time that does not depend on the size of the vector N). Alternatively, if you do this

matrix[M, N] x;


and access x[m], then there are two problems. First, we have to allocate an N-dimensional vector and copy. Second, the copy is not memory local because matrices are stored in column-major rather than row-major order.

Got it, thanks! I will switch from matrix to arrays.

Hi, may I double check, is array int T the same as int T? I am a bit confused about the variable type in Stan.

Yes, array is a newer keyword to replace that older syntax