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?

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