Posterior distribution of model variable

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?

3 Likes

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.

1 Like

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