I am trying to do integration in my Stan program.
To this end, I have the following example mode:
data {
int<lower=0> N;
vector[N] y;
real y_prior;
}
parameters {
real mu_y;
}
model {
mu_y ~ normal(y_prior, 1/N);
y ~ normal(mu_y, 1);
}
Now, at each sample iteration, I would like to find the CDF/integrate y up to some value (say 0 for simplicity).
Is there any way to do this?
Can one simply add:
generated quantities {
real y_integ;
y_integ = normal_cdf(0, mu_y, 1);
}
?
Some experimenting has made me question if this is the correct way of doing it or not.
The model:
\mu_y \sim N(y_{prior}, 1/N)
y \sim N(\mu_y, 1)
Let \Phi_y be the CDF of y include the following term in the sampling:
y_{integ} = \Phi_y(0)
The definition of normal_cdf
is,
\textrm{normal_cdf}(v \mid \mu, \sigma) = \displaystyle \int_{-\infty}^v \textrm{normal}(z \mid \mu, \sigma) \, \textrm{d}z.
Then \Phi(x) is just \textrm{normal_cdf}(x \mid 0, 1), so that we can also formulate this as
\displaystyle \textrm{normal_cdf}(v \mid \mu, \sigma) = \Phi\left(\frac{v - \mu}{\sigma} \right).
If you define that in generated quantities as you have done, then y_integ
will get a posterior distribution of CDF values based on sampled values of mu_y
. In this case, it’s a transform of the random variable mu_y
.
Ok, thank you! I was not sure this was correct but apparently, it was. Thank you for confirming and sorry if this seemed like a trivial question.