Hello,

Thank you very much for your help and your time.

I simplified the problem to two toy examples, the first one doesn’t include any generation of sub-vectors within the model, but it shows that he code works (at least sometimes as this is a toy example).

The objective is to have some IHPP of regression type with error in variables (Facotor_1). I can send R files if needed.

rm(list = ls())

library(rstan)

# Simulate some dummy data --------------------------------------------------

N0=500

N1=150

Factor1_mu_0=runif(N0,1,2)

Factor1_sd_0=rep(0.2,length(Factor1_mu_0))

IDX_1=sort(sample(1:N0,N1))

Factor1_1=Factor1_mu_0[IDX_1]

# Toy Stan model - working -----------------------------------------------------

MODEL_SPEC <- ’

functions {

real logLikelihood (

int N0,

int N1,

vector Factor1_0,

vector Factor1_1,

real beta

)

{

real logLik;

logLik=sum(Factor1_1*beta)-sum(beta*Factor1_0);

return (logLik);

}

}

data {

int<lower=0> N0;

int<lower=1> N1;

vector[N0] Factor1_mu_0;

vector[N0] Factor1_sd_0;

vector[N1] Factor1_1;

}

parameters {

real<lower=0> beta;

vector[N0] Factor1_0;

}

model {

for(i in 1:N0){

Factor1_0[i]~normal(Factor1_mu_0[i], Factor1_sd_0[i]);

}

beta~normal(0,2);

target +=logLikelihood (N0,N1,Factor1_0,Factor1_1,beta);

}

’

MODEL = stan_model(model_code = MODEL_SPEC, model_name=“stanmodel”)

MODEL_SIM = sampling(MODEL,

data = list(N0 = N0, N1 = N1, Factor1_mu_0= Factor1_mu_0,Factor1_sd_0=Factor1_sd_0, Factor1_1= Factor1_1),

iter = 5000, warmup = 1000, chains = 1, cores = 1,control = list(adapt_delta = 0.80))

MODEL_POST = extract(MODEL_SIM)

par(mfrow=c(1,2))

hist(MODEL_POST$beta,xlab= expression(beta),main=’’)

plot(MODEL_POST$beta,xlab= expression(beta),main=’’,ylab=‘trace’, type=‘l’)

and here’s the model that includes sub-vectors within the model (i excluded form the code simulation of the data):

# Toy Stan model -not working -----------------------------------------------

MODEL_SPEC <- ’

functions {

real logLikelihood (

int N0,

int N1,

vector Factor1_0,

int<lower = 1, upper = N1> IDX_1, // 1 difference

real beta

)

{

row_vector[N1] Factor1_1; // 2 difference

Factor1_1=Factor1_0[IDX_1]; // 3 difference

real logLik;

logLik=sum(Factor1_1*beta)-sum(beta*Factor1_0);

return (logLik);

}

}

data {

int<lower=0> N0;

int<lower=1> N1;

vector[N0] Factor1_mu_0;

vector[N0] Factor1_sd_0;

int<lower = 1, upper = N1> IDX_1; // 4 difference

}

parameters {

real<lower=0> beta;

vector[N0] Factor1_0;

}

model {

for(i in 1:N0){

Factor1_0[i]~normal(Factor1_mu_0[i], Factor1_sd_0[i]);

}

beta~lognormal(0,2);

target +=logLikelihood (N0,N1,Factor1_0,IDX_1,beta); // 5 difference

}

’

MODEL = stan_model(model_code = MODEL_SPEC, model_name=“stanmodel”)

MODEL_SIM = sampling(MODEL,

data = list(N0 = N0, N1 = N1, Factor1_mu_0= Factor1_mu_0,Factor1_sd_0=Factor1_sd_0, IDX_1= IDX_1),

iter = 5000, warmup = 1000, chains = 1, cores = 1,control = list(adapt_delta = 0.80))

MODEL_POST = extract(MODEL_SIM)

par(mfrow=c(1,2))

hist(MODEL_POST$beta,xlab= expression(beta),main=’’)

plot(MODEL_POST$beta,xlab= expression(beta),main=’’,ylab=‘trace’, type=‘l’)

Thank you very much!

Kind Regards

Peter