Quantile computations in Stan vs custom computation

Modeling
1

Hello,

I am fitting a linear model and I am having a small trouble with the quantiles computation. In particular I have:

transformed parameters {
  
  vector[N] mu = alpha + beta*x;

}

model {

  beta  ~ normal( 0, sigma_beta  );
  alpha ~ normal( mu_alpha, sigma_alpha );
  
  y ~ normal(mu, sigma);
}

and I am interested in computing the quantiles of mu. As far as I know it is correct to compute the quantiles of alpha and beta, and propagate them with the transformation. In other words I can do quantile_alpha and quantile_beta and get quantile_mu = quantile_beta * x + quantile_alpha. I am basing this claim on this paper http://www.gatsby.ucl.ac.uk/~snelson/gpwarp.pdf (sec 3.2).

However I am getting different results with this computation that with the quantiles that Stan returns. Can someone provide light into where this difference comes from?. I guess stan obtains mu_post using alpha_post and beta_post and then compute quantiles directly from mu_post

2

That paper only considers a one-dimensional warping function.
If either alpha or beta is known with certainty then the quantiles of the other are related to the quantiles of mu as you describe but I don’t think that holds when both are uncertain.

For example, let’s say you have x=1 so that mu = alpha + beta and medians of both alpha and beta are zero. Now you know that there’s a 25% chance that both alpha and beta are < 0 and thus mu < 0, and a 25% chance that both alpha and beta are > 0 and thus mu > 0 but you still can’t say anything about the remaining 50% case where the signs of alpha and beta disagree.
There’s no guarantee that median(alpha) + median(beta) = 0 is going to be the 50% quantile of mu, only that it is some quantile between 25% and 75%.

1 Like
3

Okay that is what I was hoping, that the result from that paper only holds for a one-dimensional warping function.