Hi, I need help with a user-built function. I am trying to run a model in which I have to estimate the parameter space to calculate an interaction surface, it involves creating a variable with my function surface_Stan(), which has some additional parameters to calculate (philength, philengthstrat, phistrat) in a linear regression.
I want the surface_Stan() to produce a vector that I can use as a predictor in a linear model, and to get the best posterior distribution for the parameters (philength, philengthstrat, phistrat) for which I also have priors. The function works well when I tested using expose_stan_functions(), but when I am writing the model I get the following message:
###############
YNTAX ERROR, MESSAGE(S) FROM PARSER:
Info: left-hand side variable (name=n) occurs on right-hand side of assignment, causing inefficient deep copy to avoid aliasing.
Info: left-hand side variable (name=n) occurs on right-hand side of assignment, causing inefficient deep copy to avoid aliasing.
error in âmodelf2b6abe01f8_3e9139b8f675d3a48a22d966bc4dff69â at line 175, column 22
173: a ~ normal( 0 , 1 );
174:
175: denP = surface_Stan(Tank_N, âŚ
----------------------------------^
176:
PARSER EXPECTED: â(â
Error in stanc(file = file, model_code = model_code, model_name = model_name, :
failed to parse Stan model â3e9139b8f675d3a48a22d966bc4dff69â due to the above error.
#############
functions {
int num_lt(real x, vector y) {// number of elements with tank x
int n = 0;
for (i in 1:rows(y))
if (y[i] == x)
n = n +1;
return n;
}
int[] find_lt(real x, vector y) { // gets the tanks ids like in the which funciton
int result[num_lt(x, y)];
int n = 1;
for (i in 1:rows(y)) {
if (y[i] == x) {
result[n] = i;
n = n + 1;
}
}
return result;
}
vector col_Sums(matrix X) {
vector[cols(X)] s ;
for (i in 1:cols(X)) s[i] = sum(col(X, i)) ;
return s ;
}
vector row_Sums(matrix X) {
vector[rows(X)] s ;
for (i in 1:rows(X)) s[i] = sum(row(X, i)) ;
return s ;
}
real[] p_size(int x, vector y, vector z) { // gets the size values per tank
real p_size[num_lt(x,y)];
int pos[num_lt(x,y)];
pos = find_lt(x,y);
for (i in 1:num_lt(x,y)){
p_size[i] = z[pos[i]];
}
return p_size;
}
real[] p_spec(int x, vector y, vector spec) { // gets the size values per tank
real p_spec[num_lt(x,y)];
int pos[num_lt(x,y)];
pos = find_lt(x,y);
for (i in 1:num_lt(x,y)){
p_spec[i] = spec[pos[i]];
}
return p_spec;
}
vector p_alpha(int x, vector y, vector spec, vector z,
real philength, real philengthstrat, real phistrat,
real equalspplength) { // gets the size values per tank
real sz[num_lt(x,y)];
real sp[num_lt(x,y)];
vector[num_lt(x,y)] alphaN;
matrix[num_lt(x,y),num_lt(x,y)] complength;
sz = p_size(x, y, z);
sp = p_spec(x, y, spec);
for ( j in 1:num_lt(x,y)){
for ( k in 1:num_lt(x,y)){
if (sp[j] == 1 && sp[k]==1) {
complength[j,k] = exp( ( philength + philengthstrat) * (sz[k] - equalspplength ) - ( philength +
philengthstrat ) * ( sz[j] - equalspplength ) );
}else if(sp[j] == 1 && sp[k]==0){
complength[j,k] = exp( -phistrat ) * exp( (philength)*( sz[k] - equalspplength ) - ( philength +
philengthstrat) *( sz[j] - equalspplength ) );
}else if(sp[j] == 0 && sp[k]==0){
complength[j,k] = exp( ( philength ) * ( sz[k] - equalspplength ) - ( philength ) * ( sz[j] - equalspplength ) );
}else if(sp[j] == 0 && sp[k] == 1){
complength[j,k] = exp( phistrat ) * exp( ( philength + philengthstrat )*( sz[k] - equalspplength ) -
( philength ) * ( sz[j] - equalspplength ) );
}
}
}
alphaN[1:num_lt(x,y)] = row_Sums(complength);
return alphaN;
}
vector surface_Stan(int Tank_N, int N, vector y, vector spec, vector z,
real philength, real philengthstrat, real phistrat,
real equalspplength){
vector[N] long_version;
int x;
for ( i in 1:Tank_N){
x = i;
if ( num_lt(x,y) > 1){
long_version[find_lt(x,y)] = p_alpha(x, y, spec, z, philength, philengthstrat, phistrat, equalspplength);
}
}
return long_version;
}
}
data{
int<lower=1> N;
vector[N] growth;
vector[N] spec;
vector[N] z;
vector[N] z2;
int Tank_N;
real equalspplength;
}
parameters{
real a;
real b_zs;
real b_g;
real b_z;
real b_z2;
real b_N;
real b_z2s;
real b_AS;
real philength;
real phistrat;
real philengthstrat;
real<lower=0> sigma;
}
model{
vector[N] gamma_3;
vector[N] gamma_2;
vector[N] gamma_1;
vector[N] mu;
vector[N] denP;
sigma ~ cauchy( 0 , 1 );
philengthstrat ~ normal( 0 , 1 );
phistrat ~ normal( 0 , 1 );
philength ~ normal( 0 , 1 );
b_AS ~ normal( 0 , 1 );
b_z2s ~ normal( 0 , 1 );
b_zs ~ normal( 0 , 1 );
b_N ~ normal( 0 , 1 );
b_z2 ~ normal( 0 , 1 );
b_z ~ normal( 0 , 1 );
b_g ~ normal( 0 , 1 );
a ~ normal( 0 , 1 );
denP = surface_Stan(Tank_N = Tank_N, N = N, y = tank, spec = spec, z= z, philength, philengthstrat, phistrat, equalspplength);
for ( i in 1:N ) {
gamma_3[i] = b_N + b_AS * spec[i];
}
for ( i in 1:N ) {
gamma_2[i] = b_z2 + b_z2s * spec[i];
}
for ( i in 1:N ) {
gamma_1[i] = b_z + b_zs * spec[i];
}
for ( i in 1:N ) {
mu[i] = a + b_g * spec[i] + gamma_1[i] * z[i] + gamma_2[i] * z2[i] + gamma_2[i] * denP[i];
}
growth ~ normal( mu , sigma );
}
generated quantities{
vector[N] gamma_3;
vector[N] gamma_2;
vector[N] gamma_1;
vector[N] mu;
real dev;
dev = 0;
for ( i in 1:N ) {
gamma_3[i] = b_N + b_AS * spec[i];
}
for ( i in 1:N ) {
gamma_2[i] = b_z2 + b_z2s * spec[i];
}
for ( i in 1:N ) {
gamma_1[i] = b_z + b_zs * spec[i];
}
for ( i in 1:N ) {
mu[i] = a + b_g * spec[i] + gamma_1[i] * z[i] + gamma_2[i] * z2[i] + gamma_2[i] * denP[i];
}
dev = dev + (-2)*normal_lpdf( growth | mu , sigma );
}
``
data_list <- list(
N <- dim(d)[1],
tank <- d$tank,
spec <- d$spec,
z <- d$z,
z2 <- d$z2,
growth <- d$growth,
Tank_N <- 68,
equalspplength= 18
)
Mod.1 <- stan(model_code = Stan.mod.1, data = data_list, init = 1e4, chains = 1)
####################
Thank you for your help, please find the data in this link :
https://docs.google.com/spreadsheets/d/1VTGJwD91SwPeS5p6KaoT1CSJQMd33dMIhyF7qTSCeL0/edit?usp=sharing
Jaime