I have read through other questions on the site and I feel that none provide a great answer to this question. Simply put - could anyone explain, common approaches for generating predictions on new unseen data with a Bayesian model.

For example, suppose you have a simple linear regression: $$ y = \alpha + \beta * x $$

Which under a probabilistic framework is given by:

$$ y \sim N(\mu = \alpha + \beta *x, \sigma = \epsilon) $$

We set prior distributions for the parameters $ \alpha, \beta$ and $\epsilon$ and with our training data $X_{train}$ we return samples from posterior to obtain posterior distributions across these parameters.

My question is, given new data $X_{new}$, how are predictions then made on that data?

For example, would you:

  • For each new data point in $X_{new}$, take the posterior samples from the 'training' process, and calculate $Y_{new}$ using combinations of coefficients from the posterior samples (for a given data point in $X_{new}$) - to generate a vector of predictions for each data point in $X_{new}$?


You have draws from posterior, so for a new $x_i$ the corresponding prediction for $\mu$ is $\dfrac{1}{N} \sum _i \alpha_i + \beta_i x%$. If you're familliar with rstanarm, this is what posterior_linpred does. Since the predictor is linear, this should be equal to $E(\alpha) + E(\beta)x$, where $E$ is the expectation. Here is an example using Stan.

Let's first set up the model and fit it.


model_code = '
  int N;
  vector[N] x;
  vector[N] y;
  real alpha;
  real beta;
  real<lower=0> sigma;
  alpha ~ normal(0,1);
  beta ~ normal(0,1);
  sigma ~ cauchy(0,1);
  y ~ normal(alpha + beta*x, sigma);
generated quantities{
  real yppc[N] = normal_rng(alpha + beta*x, sigma);
fl = write_stan_file(model_code)
model = cmdstan_model(fl)

N = 10
x = rnorm(N)
y = 2*x + 1 + rnorm(N, 0, 0.45)
model_data = list(N=N, x=x, y=y)
fit = model$sample(model_data, parallel_chains=4)

Now, extract the draws

beta = as.numeric(fit$draws('beta'))
alpha = as.numeric(fit$draws('alpha'))
sigma = as.numeric(fit$draws('sigma'))

Let's say I want to make a prediction for $x=3$. I should compute $\alpha_i + \beta_i \cdot 3$ for all my posterior $\alpha_i, \beta_i$. In R...

prediction = alpha + beta*3

enter image description here

Note that we get a distribution of possible values. Our prediction in this case is the mean of this distribution


Which is close to the true value of 7.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.