Bayesian Prediction Simple Explanation 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}$?

 A: Correct.
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.
library(tidyverse)
library(tidybayes)
library(cmdstanr)

model_code = '
data{
  int N;
  vector[N] x;
  vector[N] y;
}
parameters{
  real alpha;
  real beta;
  real<lower=0> sigma;
}
model{
  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
hist(prediction)


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

Which is close to the true value of 7.
