I’m not sure what you’re suggesting.

you can’t reference variables declared in a later program block in an earlier program block, so no, you can’t set things up that way.

I’m not sure what you’re suggesting.

you can’t reference variables declared in a later program block in an earlier program block, so no, you can’t set things up that way.

In my model the `farenheit_reported`

and `celcius_reported`

are observed variables, i.e., data. However, as you have recommended their declaration goes in the `generated quantities`

block.

My question is whether I can still pass in the observed values as part of the `data`

argument in the `stan()`

call. And if I do that will my posterior be conditioned on those variables?

Apologies if I am not being very clear and thanks much for your help.

your posterior is determined by the statements in the model block.

to deal with truncated data, you should check the Stan manual.

if your `_reported`

information is coming in as data, then it should be declared in the `data`

block.

the `generated quantities`

block is run *after* the posterior has been updated. its purpose is to computed quantities of interest based on the posterior.

again, please find examples in the Stan manual or in Stan’s example models on how work with truncated data. I’m sorry I can’t answer this myself very well, as I’m a programmer with very little modeling experience.

Thanks for your help, in any case. The data is not truncated but clipped, i.e., the values outside the range are not thrown away but observed as the edge values.

I don’t quite understand how to incorporate observed data that are transformations of other variables. The few different things I tried (putting it in the model, in transformed params etc.) all gave the same error. I couldn’t find the answers in the Stan manual either.

Do you recommend posting the modeling question in the forum again as another topic?

Those are called “censored”. The terminology is terribly confusing. I go over it in the manual in the chapter “Truncated or Censored Data”. For censoring, you basically just code up a reading of min using the lcdf or the max using the lccdf.

What are you trying to do? You can define transformed data variables that are transforms of the actual data you read in. Or you can define transformed parameters that combine data and parameters (those can also be local variables in the model if you don’t want to save them).

Thanks for the tip. I will read the Stan manual to see if I can follow the examples about using censoring for my problem.

I have a more general question about modeling processes of the following type. Assume

```
z1 ~ Normal(mu1, sigma1)
z2 ~ Normal(mu2, sigma2)
x = f(z1, z2)
```

where `f`

is a deterministic function.

If I observe `x`

and am trying to infer `E[z1, z2 | x]`

, how do I go about modeling this in Stan?

Looks like I’ll have to put `x`

in the data section, but when I try doing that but declare `x = f(z1, z2)`

in the model section I get the “cannot assign variable outside declaration block” error.

This is partly a Stan question and partly a stats question.

For Stan, you need to define any derived quantities that depend on parameters in either the transformed parameters block or as a local variable in the model block (or in generated quantities if the log density doesn’t depend on the output).

The stats problem is more subtle and is going to depend on what `f`

is and other information you have in the model. Typically, you need a Jacobian adjustment when changing variables this way and trying to put a distribution on the transform. Here, you need to be able to map from `x`

back to `z1`

and `z2`

. The manual explains how this works in the chapter on transformed parameters and also the entire chapter on how the internal transform works.

Bob,

Since `f`

is not invertible, how would I specify a map from `x`

to `z1`

and `z2`

? For a simple case where `x = z1^2 + z2^2`

the log-likelihood of x is difficult to derive from the fact that `z1`

and `z2`

are normally distributed. How does Stan go about performing inference of `E[z1|x]`

? Could you suggest a reference which explains how the log-likelihood calculation works in the cases where Bayesian models have non-bijective transformations of parameters?

I couldn’t find the string “internal transform” in the Stan manual. Do you perhaps mean “inverse transform”?

chapter “Transformations of Constrained Variables”

try rewriting the model in terms of the data-generating process directly, instead of using the function `f`

to recover multiple parameters based on a single number.

Thanks for that.

Unfortunately in my case the data generation process is naturally specified as how I have in the example. My observations are outputs of a deterministic function of latent variables, and I am interested in recovering the latent variables.

The fundamental problem is degrees of freedom. If you have multiple latent variables that are run through a many-to-one function, and then you try to fit with the result, the estimates for the latent variables are not determined. As a simple example, in `a + b = 5`

, the values of `a`

and `b`

are not determined. Same thing in a statistical model if you try to model `(a + b)`

and then use that to recover `a`

and `b`

.

If you don’t have independent information on `a`

and `b`

there’s just nothing you can do. In the statistical setting, you can add priors, but then the scale of `a`

and `b`

is fully determined by the prior.

This is the fundamental reason you can’t compute the Jacobian—it needs to preserve volume (or area in the 2D case). If you go from 2 dimensions (area) down to 1 dimension (line), you no longer have any area. So there’s no Jacobian that can compensate for the change in area.

1 Like

I am precisely in the setting where I have priors. In your example assume

```
a ~ Normal(1, sd = 1)
b ~ Normal(0, sd = 2)
x = a + b
```

Then a, b and x are jointly Gaussian with mean vector `(1, 0, 1)`

and covariance matrix

1 0 1

0 4 4

1 4 5

Even though this matrix is singular, I can still make pairwise inferences. In particular I can calculate `E[a | x = 5] = 1/5 * (5 - 1) + 1= 1.8`

, which is not trivially obvious from the model and the observation that `x= 5`

. (In fact, `a|x=5 ~ Normal(1.8, sqrt(.8))`

)

The above calculation seems straightforward and completely in keeping with the principles of Bayesian reasoning (I hope there are no errors in my calculations above).

What I don’t understand it is what about automated inference in Stan (and other PPLs) makes this hard.

If a, b, and x were jointly Gaussian, then they would have support in all of R^3. But that’s not the case—they only have support where x = a + b. There’s no way to model this in a proper multivariate normal because the covariance matrix isn’t positive definite. You can see the rank degeneracy of your matrix just by inspection of the columsn (the third of which is the sum of the first two). That means you can’t invert it, which means you can’t evaluate the quadratic form inside the definition of the density.

The more general problem this induces is of identifiabiltiy. The data doesn’t determine `a`

or `b`

, just their sum, and there are infinitely many possibilities producing that sum.

Agreed. I did allude to the fact that the covariance matrix is singular in my post. However, as I showed the conditional density of `a|x = 5`

is definitely *not* degenerate.

Even though jointly a, b and x are on a plane, i.e., from a degenerate distribution, the priors on `a`

and `b`

do render the posterior of `a`

and `b`

given `x = 5`

a proper distribution.

As Gelman says in this post http://andrewgelman.com/2014/02/12/think-identifiability-bayesian-inference/ “a Bayesian model is identified if the posterior distribution is proper”

OK, then I think we’re on the same page. The likelihood isn’t classically identified, but with a prior, you get identification. The prior gives you identification even without data. Andrew has several notions of identification just in his book with Jennifer—you need to generalize to Bayes and to practice where you actually care about particular data sets.

To be clearer about what Stan’s trying to do, it unconstrains the parameters to have support on R^N. It assumes that it has support over all that, not just over a lower-dimensional submanifold. Stna will let you build the above model

```
a ~ foo(...);
b ~ bar(...);
x = a + b;
x ~ baz(...);
```

The real problem is that in this simple case, `a`

and `b`

are only identified through their priors.

Stan tries to unconstrain the parameters to have support on R^N. We don’t allow hard relations that produce submanifolds as there’s no way to produce a volume-preserving Jacobian. So we could parameterize this in terms of `a`

and `b`

as above and `x`

could be defined as a transformed parameter.

I see. Thanks. I assume `baz()`

in your example is the marginal distribution of `x`

, derived from the fact that `x = a+b`

.

That might work for my problem, where I have an observed quantity that is a deterministic non-linear function of some latent parameters (for which have priors). I will just have to figure out how to derive the marginal distribution of my observed quantity.

ps. I really appreciate your (and Mitzi’s) thoughtful responses to my questions.

That’s how it usually comes up. Stan lets you just define `x = a + b`

as a local variable or transformed parameter if it’s easier to put distributions on `a`

and `b.

Can I just clarify that in the situation where we define `x = a + b`

in the transformed parameters block and then put a likelihood (I’m not sure if that’s what it should be called) on x in the parameters block as `x ~ baz()`

, we do not have to add a Jacobian or similar adjustment to the model likelihood?

In my case I would like to constrain a derived parameter `x = a * b / (c * d)`

. In this case, is what is described in this thread still relevant?

Yes. If you have `a`

,`b`

,`c`

, and `d`

as `parameters`

and `x = f(a,b,c,d) ~ baz()`

in `model`

then you’ll face exactly the same problem, namely there’s no Jacobian adjustment for f and the likelihood is not sufficient to identify the model.

1 Like

Hi there,

I am pretty new to R/RStan, so apologies in advance if this doesn’t make much sense. I have uploaded the RStan file I have been using to this post. When I run the following R code:

library(readr)

library(“rstan”)

options(mc.cores=parallel::detectCores())

rstan_options(auto_write = TRUE)

setwd(“c:/Users/jlevy/Dropbox (Sydney Uni)/Twins reference point/Analysis_XW/twin data/BHM”)

library(readxl)

sample ← read_csv(“/Users/jlevy/Dropbox (Sydney Uni)/Twins reference point/Analysis_XW/twin data/BHM/for_rstan_JL.csv”)

N=nrow(sample)

nsubj=max(sample[‘sid’])

subjs=unlist(sample[‘sid’])

choices=unlist(sample[‘chosea’])

x1a=unlist(sample[‘x1a’])

x2a=unlist(sample[‘x2a’])

x3a=unlist(sample[‘x3a’])

x4a=unlist(sample[‘x4a’])

x1b=unlist(sample[‘x1b’])

x2b=unlist(sample[‘x2b’])

x3b=unlist(sample[‘x3b’])

x4b=unlist(sample[‘x4b’])

p1a=unlist(sample[‘p1a’])

p2a=unlist(sample[‘p2a’])

p3a=unlist(sample[‘p3a’])

p4a=unlist(sample[‘p4a’])

p1b=unlist(sample[‘p1b’])

p2b=unlist(sample[‘p2b’])

p3b=unlist(sample[‘p3b’])

p4b=unlist(sample[‘p4b’])

sid=unlist(sample[‘sid’])

sq=unlist(sample[‘sq’])

eva=unlist(sample[‘eva’])

evb=unlist(sample[‘evb’])

maxmin=unlist(sample[‘maxmin’])

minmax=unlist(sample[‘minmax’])

xatmaxp=unlist(sample[‘xatmaxp’])

fit<- stan(file=‘bhm_5rp_JL.stan’,

data=list(N=N,nsubj=nsubj,subjs=subjs,choices=choices,x1a=x1a,x2a=x2a,x3a=x3a,x4a=x4a,

x1b=x1b,x2b=x2b,x3b=x3b,x4b=x4b,p1a=p1a,p2a=p2a,p3a=p3a,p4a=p4a,p1b=p1b,p2b=p2b,p3b=p3b,p4b=p4b,sid=sid,sq=sq,eva=eva,evb=evb,maxmin=maxmin,minmax=minmax,xatmaxp=xatmaxp),

warmup=2500,

chains=4,

iter=10000,

init_r=1,

verbose=FALSE)

print(fit)

fit_ss ← extract(fit,permuted=TRUE)

I get the following error message:

**Error in stanc(file = file, model_code = model_code, model_name = model_name, : **
** parser failed badly**

Can someone please advise?

Thanks in advance!

bhm_5rp_JL.stan (18.3 KB)