Regression in Baysian settings Assume we have the posterior distribution of this linear regression model $y = w^Tx$, $P(w | D,\theta)$, where $D = \{(x_i,y_i)\}_{i \in \{1,\dots,n\}}, n $ is the number of data instances, $\theta$ is the set contains $w$ and all hyperparmaters. For a testing point $\hat{x}$, we calculated the posterior distribution as following:
$$P(w|D,\theta) = \int_\theta P(\hat{x} | \theta) P(\theta|D)d\theta$$.
Since we have a continous possible values of $w$, how to "discretize" the values of $w$?
Assume we already did that, and we have several values of $w$, $\{w_1,\dots,w_k\}$, with $P(w = w_j | D,\theta), j \in \{1,\dots,k\}$, then the prediction based on these different values of $w_j$ must be associated with their posterior values. For instance, I'd trust the prediction $y_1$ (i.e., based on $w_1$) more than $y_3$  as $P(w = w_1 | D,\theta) > P(w = w_3 | D,\theta)$. How to reflect such confidence to the predicted $y$?

A common approach is to to approximate the posterior (as oppose of calculating the full posterior). We usually model the posterior with "easy to work with distribution" for example distribution that can be described by finite number of parameters.
A common solution is to model the posterior as a normal distribution, and thus for each parameter we only need to find $\mu, \sigma$. In this case, the prior is also modeled as the standard normal distribution, and the way to learn the (approximated) posterior is via variational inference.
Once you learn the posterior you can do prediction by computing the integral of the posterior distribution i.e. weighted average according to the posterior weights. Practically, computing the integral is impossible (as there are infinite values) and thus the average is done by sampling the posterior (in this setup sampling each parameter according its normal distribution.
