Dear Forum,
What is the most efficient way to retrieve the posterior distribution of a model variable, i.e. a compound variable that is declared and calculated in the model section.
Thanks,
Christian Damgaard
Dear Forum,
What is the most efficient way to retrieve the posterior distribution of a model variable, i.e. a compound variable that is declared and calculated in the model section.
Thanks,
Christian Damgaard
Hmm, I think what you want is to move the declaration/definition of the variable to the transformed parameters
block:
Any variable declared as a transformed parameter is part of the output produced for draws.
Contrast this with the model
block:
The variables in the model block are local variables and are not written as part of the output.
Does this help?
Thanks. It works.
Christian
It works, but we have to duplicate most of the model into the generated quantities, which is awkward. Why is it not possible to declare a global variable / parameter that record the content of the local variable in the model section?
Christian
@Funko_Unko didnât suggest moving the variable to generated quantities
, but rather to transformed parameters
. The transformed parameters can be reused in the model block.
Sorry. It actually did not work when we declared the latent variable in âtransformed parametersâ, then the model section compiled with an error. We had to define the latent variable in the model section and then calculate the latent variable again in âgenerated quantitiesâ. Which is awkward.
If you post the code that didnât work and the code that did, I bet somebody here can spot the bug in the code that didnât work.
The same code should work in any block. Iâm guessing it didnât work because you moved only the declaration to transformed parameters
but tried to assign the value in model
block. You need to move all the calculations to the transformed parameters
block.
Here is the code that works
data {
int<lower=0> n; // number of data rows in Measurements
int<lower=0> nNP; // number of data rows in LoadNP
vector[n] chl; // chlorofyl measurements
matrix[nNP,2] nitphos; // N and P loadings
//Fixed effects (intercept, BV, Sali, Temp, âŚ)
int<lower=0> p; // Colums in design matrix
matrix [n,p] X; // Design matrix
//Hierarchical structure
int<lower=0> j; // number of stations
int<lower=0> k; // number of years
int<lower=1,upper=j> stationID[n]; // vector of integer that assign each data row in Measurements to a station
int<lower=1,upper=k> yearID[n]; // vector of integer that assign each data row in Measurements to a year
int<lower=1,upper=j> NPstationID[nNP]; // vector of integer that assign each data row in LoadNP to a station
int<lower=1,upper=k> NPyearID[nNP]; // vector of integer that assign each data row in LoadNP to a year
}
parameters {
vector[j] a; // Parameter for effects of N and P
vector[j] b; // Parameter for effects of N and P
vector[j] c; // Parameter for effects of N and P
vector[p] beta; // Parameter for effects of factors in X
cov_matrix[2] Sigma; // Covariance matrix for N and P loadings
real<lower=0> chlscale; // Scale parameter for chlorofyl measurements
//Random effect
vector[j] etarandom;
real<lower=0> sdrandom;
//Process error
vector[n] etaproc;
real<lower=0> sdproc;
//Latent variables
matrix[j*k,2] NPtrue;
}
transformed parameters {
// vector[n] chltrue; If chltrue is declared here, then the model section does not compile
}
model {
int jn;
int kn;
matrix[nNP,2] NPtrueforallNP;
vector[n] Ntrue;
vector[n] Ptrue;
vector[n] NandP;
vector[n] mu;
vector[n] randrandom;
vector[n] randproc;
vector[n] chltrue;
vector[n] alfagamma;
vector[n] betagamma;
//Priors
a ~ normal(0, 10000);
b ~ normal(0, 10000);
c ~ normal(0, 10000);
beta ~ normal(0, 10000);
Sigma ~ inv_wishart(2, diag_matrix(rep_vector(0.1,2)));
chlscale ~ exponential(0.1);
etarandom ~ normal(0, 1);
sdrandom ~ exponential(0.1);
etaproc ~ normal(0, 1);
sdproc ~ exponential(0.1);
col(NPtrue, 1) ~ normal(5, 100);
col(NPtrue, 2) ~ normal(5, 100);
//Models
for (iNP in 1:nNP){
jn=NPstationID[iNP];
kn=NPyearID[iNP];
NPtrueforallNP[iNP]=NPtrue[(jn-1)*kn+kn];
}
for (i in 1:n){
jn=stationID[i];
kn=yearID[i];
Ntrue[i]=NPtrue[(jn-1)*kn+kn,1];
Ptrue[i]=NPtrue[(jn-1)*kn+kn,2];
NandP[i]=a[jn]*Ntrue[i]+b[jn]Ptrue[i]+c[jn]Ntrue[i]Ptrue[i]; // a(Ntrue[i]+bPtrue[i])+c(Ntrue[i]/(Ntrue[i]+Ptrue[i]))
}
mu â X*beta;
for (i in 1:n)
randrandom[i]=sdrandom*etarandom[stationID[i]];
randproc = sdproc*etaproc;
chltrue = NandP + mu + randrandom + randproc;
for (i in 1:n){
alfagamma[i] = pow(chltrue[i],2.0)/chlscale;
betagamma[i] = chltrue[i]/chlscale;
}
//Likelihood functions
for (iNP in 1:nNP)
nitphos[iNP] ~ multi_normal(NPtrueforallNP[iNP], Sigma);
chl ~ gamma(alfagamma,betagamma); //alfa = mu^2/scale, beta = mu/scale
}
@nhuurre nailed it. Move this:
chltrue = NandP + mu + randrandom + randproc;
To the transformed parameters block. Youâll also need to move all of the parameters that underly the computation of chltrue
up to the transformed parameters block, and the code to compute their values.
I had to put some more code into the transformed parameters block. But now it works as I was hoping for.
Thanks to you all,
Christian