I am trying to use the posterior predictive distribution of my model to estimate the probability that a future observation is greater than 15 for each of (1) x_new and (2) y_new.
My Stan code is as follows:
data {
int<lower=1> n;
vector[n] x;
vector[n] y;
}
parameters {
real x_mu;
real y_mu;
real<lower = 0> sigma;
}
transformed parameters {
real prob_positive_diff;
prob_positive_diff = y_mu >= x_mu;
}
model {
x ~ normal(x_mu, sigma);
y ~ normal(y_mu, sigma);
sigma ~ normal(0, 10);
}
generated quantities{
real x_new;
real y_new;
x_new = normal_rng(x_mu, sigma);
y_new = normal_rng(y_mu, sigma);
}
My understanding is that what I’ve done in the generated quantities block is sample from the posterior predictive distribution, right? I’m not 100% sure I’ve done this correctly, so I’d appreciate it if people could please provide feedback.
What I’m unsure of now is how to estimate the probability that each of (1) x_new and (2) y_new are greater than 15? I tried something like the following:
I mean that we want the two separate probabilities of x_new > 15 and y_new > 15, where x_new and y_new are samples of the posterior predictive distribution. It would be similar to prob_positive_diff = y_mu >= x_mu; in my code, except there would be 2 of them, one for each of x_new > 15 and y_new > 15.
I don’t think the code you posted samples from the posterior predictive distribution? Or perhaps I am misunderstanding something? I was under the impression that we get the posterior predictive distribution from the code I posted in the generated quantities block? I used the following examples to guide me with regards to the posterior predictive distribution:
And it seems here that you would be multiplying the probabilities:
p = exp(log1m(normal_cdf(15, mu_x, sigma)) + log1m(normal_cdf(15, mu_y, sigma)));
?
It is possible that I am thinking about this implementation of the posterior predictive distribution incorrectly, which is why I’m struggling to find out how to calculate the aforementioned probabilities that I am seeking.
Sorry, I misread your post and thought that you were looking for the probability that BOTH y_{new}>15 AND x_{new}>15, so P(y_{new}>15)P(x_{new}>15).
You probably meant this:
stancode_alt:
data{
int N;
vector[N] x;
vector[N] y;
int N_rep;
}
parameters{
real mu_x;
real mu_y;
real<lower=0> sigma;
}
transformed parameters{
real<lower=0,upper=1> p_x = exp(log1m(normal_cdf(15, mu_x, sigma)));
real<lower=0,upper=1> p_y = exp(log1m(normal_cdf(15, mu_y, sigma)));
}
model {
x ~ normal(mu_x, sigma);
y ~ normal(mu_y, sigma);
sigma ~ normal(0, 10);
}
generated quantities{
real p_x_new;
real p_y_new;
{
vector[N_rep] x_new_ind;
vector[N_rep] y_new_ind;
for(n in 1:N_rep){
x_new_ind[n] = normal_rng(mu_x, sigma) > 15 ? 1.0 : 0.0;
y_new_ind[n] = normal_rng(mu_y, sigma) > 15 ? 1.0 : 0.0;
}
p_x_new = sum(x_new_ind)/N_rep;
p_y_new = sum(y_new_ind)/N_rep;
}
}
samples_alt <- stan(model_code = stan_code_alt, data = list(N=N,x=x,y=y,N_rep=1e5))
Inference for Stan model: 9f1eca9ece6b6dc4ba656cac787e03b5.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu_x 2.79 0.02 0.96 0.94 2.15 2.79 3.43 4.68 3366 1
mu_y 5.95 0.02 0.94 4.15 5.30 5.93 6.59 7.75 3886 1
sigma 9.42 0.01 0.47 8.54 9.09 9.41 9.73 10.37 4000 1
p_x 0.10 0.00 0.02 0.06 0.08 0.10 0.11 0.14 3594 1
p_y 0.17 0.00 0.03 0.12 0.15 0.17 0.19 0.23 3890 1
p_x_new 0.10 0.00 0.02 0.06 0.08 0.10 0.11 0.14 3581 1
p_y_new 0.17 0.00 0.03 0.12 0.15 0.17 0.19 0.23 3903 1
lp__ -546.62 0.03 1.23 -549.75 -547.18 -546.31 -545.72 -545.21 2053 1
Samples were drawn using NUTS(diag_e) at Fri May 18 16:20:20 2018.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
You can play around with the values for N_rep. The higher N_rep the closer p_x_new will be to p_x, because the posterior predictive is equal to the posterior here. Right?
The higher N_rep the closer p_x_new will be to p_x, because the posterior predictive is equal to the posterior here. Right?
Yes. Assuming my understanding of the theory is correct, then, if the model (and, therefore, the posterior distribution) is accurate, then the posterior predictive should (roughly) estimate the posterior distribution itself. In this (mock/example) case, since our observations are realisations from a normal distribution, and since our simulated replicated data are also realisations from a normal distribution, we would expect them to be approximately the same. Does this sound right?
So, given this definition, it seems like normal_cdf(15, mu_y, sigma) would be the probability that it is less than or equal to 15? Or am I misunderstanding this?
This was my initial thought. Would there be any difference between this and what @Max_Mantei did in his post? In other words, both methods would be correct, right? His solution seems like it would be more elegant.
Note that log1m(x) is \log(1-x), so that exp(log1m(normal_cdf(15, mu_x, sigma))) effectively becomes 1-\int_0^{15}\text{Normal}(\mu_x,\sigma), which is probably also ok to put in the model as 1 - normal_cdf(15, mu_x, sigma). (I used the exp(log(...)), because I initially thought [wrongly] that you wanted to multiply the two probability and it’s usually better to add them on a log-scale.)
Well, actually you pointed at this (using the cdf) in your first post. I think it depends on the specifics of the application. In this case using a transformed parameter made sense. But If you have new data, that doesn’t enter your model, then perhaps the R or generated quantities approach is better, where better means ease of implementation.
Also, it depends if you want to have a distribution over possible p_x or just a number: mean(x_new > 15) will only give you one number, but sometimes that’s all you need and then this would be the easiest way to get what you want.
Yes, I’m only interested in the probability, not the distribution. In that case, I’ll leave it as is, since both implementations would get me the same information.
