What Does Stan Calculate in case of Inconsistent Modeling Assumptions

Modeling
1

I am not sure what Stan calculates in the following case:

  m <- rstan::stan_model(model_code = '
                            data{real x;}
                            parameters {real <lower=0>t;}
                            model {x ~ normal(t,1);
                         x ~ normal(t^2,1);
                         }')
  f <- rstan::sampling(m, data=list(x=1), iter = 1100,chains=1)
  
  
  stan_trace(f)
   stan_dens(f)

where, the data x is assumed that it is distributed by different two models, namely, one is

X \sim Normal(\theta,1),
and the other is

X \sim Normal(\theta^2,1).

The model converged in R hat criterion.

2

The sampling statements are translated to increments of the log density under the hood. Check out the syntax here: 7.4 Sampling Statements | Stan Reference Manual

So it’s perfectly valid Stan code to do things like:

x ~ normal(t, 1);
x ~ normal(t^2, 1);

Like you said though, it’s confusing what that means in terms of statistics.

That says that x comes from a distribution of which the density is proportional to the product of two normal densities (where the parameter is connected somehow). Presumably it’s some sorta distribution.

Rhat only checks that estimates within a chain (the front vs. the back) are comparable and between chain estimates are comparable as well. It isn’t made to tell you if your model is good or bad.

1 Like
3

Probably the way to to do that is make predictions. To make predictions with this model though you’d need to figure out what the product of those two distributions are so you can draw random numbers from it. Maybe it’s another normal with a weird mean? But I do not know.

4

Thank you for reply.

This problem has been bothering me a while. Now I am relieved. Thank you.

I think that it looks tough. :-'D

The convolution of two Gaussians is also Gaussian, namely N_{m+m',v+v'}(x)=\int_{-\infty}^{\infty} N_{m,v}(y)N_{m',v'}(x-y)dy, with simple algebraic relations between their mean and variance. But in the above code, it is not a convolution but a product of Gaussians, so I am not sure but I guess the product of two Gaussian is not so.

5

I found this answer that the product of two Gaussian random variables is distributed, in general, as a linear combination of two Chi-square random variables: https://math.stackexchange.com/a/397716

1 Like
6

Thank you! I did not know it.

7

This is good info. In this case, however, one does not have two random variables. That would be the case if you had X \sim \text{normal}(t, 1) and Y \sim \text{normal}(t^2, 1) and you were interested in Z = XY. Here one is assuming a density on X such that f_X(x) = \phi(x; t, 1)\phi(x; t^2, 1), where \phi(\cdot; \mu, \sigma) is the normal density with parameters \mu and \sigma. In @Jean_Billie’s post above f_X is a weird density. If you were to write:

target += 0.5 * normal_lpdf(x | t, 1);
target += 0.5 * normal_lpdf(x | square(t), 1);

Then X would be normally distributed with parameters \mu^\ast = \frac{t + t^2}{2} and \sigma^\ast = 1. This is called geometric pooling.

5 Likes
8

Bromiley is everybody’s friend here: http://www.tina-vision.net/docs/memos/2003-003.pdf

4 Likes