I want to add loglikelihood to generated quantities to my stan code :
functions {
real E_Riess(real zz, real om,real w,real w_a) {
return 1/ sqrt(om*(1+zz)^3
+(1-om)*exp(-3*w_a/(1.0+zz))*(1+zz)^(3*(1+w+w_a)));
}
real loglikelihood_lpdf(vector mu, vector muth, matrix invCov) {
real logL;
logL = ((mu-muth)')*(invCov*(mu-muth));
logL = logL*-0.5;
return logL;
}
}
data {
int<lower=1> nobs;
vector[nobs] mu;
vector[nobs] z;
cov_matrix[nobs] invCov;
}
transformed data {
}
parameters{
real<lower=0, upper=1.0> om;
real<lower=0> sigint;
real<lower=-2, upper=.1> w;
real<lower=-2, upper=.1> w_a;
}
transformed parameters{
vector[nobs] mu_th;
for (i in 1:nobs) {
mu_th[i] = E_Riess(z[i],om,w,w_a);
}
}
model {
target += loglikelihood_lpdf(mu| mu_th, invCov);
}
generated quantities{
vector[nobs] log_lik;
for (i in 1:nobs){
log_lik[i] =loglikelihood_lpdf(mu[i]| mu_th[i],invCov);
}
}
When I run this code I get this error:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
No matches for:
loglikelihood_lpdf(real, real, matrix)
Available argument signatures for loglikelihood_lpdf:
loglikelihood_lpdf(vector, vector, matrix)
error in 'modelf6c62144cd7_0466e8e40e3cc2338230babdeb5630ac' at line 43, column 60
-------------------------------------------------
41:
42: for (i in 1:nobs){
43: log_lik[i] =loglikelihood_lpdf(mu[i]| mu_th[i],invCov);
^
44: }
It seems like you are defining loglikelihood_lpdf inside another function, so it wouldn’t be available outside of that scope. Try specifying it as a standalone function, which I think is the best practice for writing functions anyway (although in this case the E_Riess function is called only once so it wouldn’t make a practical difference).
Thanks for your reply. But if I used stan function multi_normal_lpdf instead of loglikelihood_lpdf I get the same error. What I don’t understand is why both arguments assign real even though I declared both as vectors.
In generated quantities you are running your function for each i separately. You defined the lpdf to ingest vectors, but line 43 passes reals rather than vectors. You need to write a different function that ingests reals and is not vectorised in your case.
I suggest something like;
real loglikelihood_lpdf(int i, real mu, real muth, matrix sigma) {
return normal_lpdf(mu|muth,sigma[i,i]);
}
as far as I can tell, you want the marginal likelihood. but i suggest that this makes little sense if you are using a multivariate normal. In that case the log_lik should be defined as a scalar because you only get a loglikelihood at each point in your n dimensional space (nobs number of dimensions). So only one loglikelihood for each nobs dimensional vector.
Yes. As the message says, it is about the signature. It requires type vector for the first two arguments. You can probably fix it by using mixed operation functions like to_vector.
potentially, you are asking for an evaluation of the lpdf at a point in parameter space that assumes you are at the mean of the multidimensional normal for all observations bar one i.e vector_x[j != i] = mu_th[j], vector_x[i] = mu[i]. then log_lik = loglikelihood_lpdf(vector_x | mu_th,invCov);
It is just a little unclear what you are trying to achieve in the generated quantities section.
Thanks for your reply. I tried to use your function but unfortunately, I get an error. without generated quantities block my code run correctly but for model checking and leave one out cross-validation approach, I need to calculate loglikelihood. Here is my code and data:
library(rstan)
Cov_inv = read.table(invcovRiess.dat'))
zmu_data = read.table(Riessinv.dat'))
data1 = list(nobs = 6, z = zmu_data[,1], mu = zmu_data[,2],invCov = Cov_inv)
model_code1= "
functions {
real E_Riess(real zz, real om,real w,real w_a) {
return 1/ sqrt(om*(1+zz)^3
+(1-om)*exp(-3*w_a/(1.0+zz))*(1+zz)^(3*(1+w+w_a)));
}
real loglikelihood_lpdf(vector mu, vector muth, matrix invCov) {
real logL;
logL = ((mu-muth)')*(invCov*(mu-muth));
logL = logL*-0.5;
return logL;
}
}
data {
int<lower=1> nobs;
vector[nobs] mu;
vector[nobs] z;
cov_matrix[nobs] invCov;
cov_matrix[nobs] Cov;
}
transformed data {
}
parameters{
real<lower=0, upper=1.0> om;
real<lower=0> sigint;
real<lower=-2, upper=.1> w;
real<lower=-2, upper=.1> w_a;
}
transformed parameters{
vector[nobs] mu_th;
for (i in 1:nobs) {
mu_th[i] = E_Riess(z[i],om,w,w_a);
}
}
model {
target += loglikelihood_lpdf(mu| mu_th, invCov);
}
generated quantities{
}
"
fit <- stan(model_code = model_code1,
data = data1,
seed = 42,
chains = 1,
iter =1000,
cores= 1,
warmup=500,control = list(adapt_delta = 0.99999)
)
print(fit,pars=c("om","w","w_a"),intervals=c(0.025, 0.975), digits=3)
my function was not designed to work. it was to provide guidance to you regarding what you are trying to achieve in generated quantities.
Now that you have removed generated quantities, there is no issue with your code. Your loglikelihood_lpdf function calculates the loglikelihood as you desired and it is stored in the model results as you need.
the loglikelihood is returned for each MCMC step as part of the final results and you can extract it several ways e.g. the extract() function.
So you don’t need the LL to be calculated in the generated quantities block as far as I can tell.