TL;DR to "how does maximizing ELBO lead to a good/correct posterior predictive distribution 𝑝(𝑦|𝑥,𝐷)?"
The initial formulation for making new predictions is correct, essentially you have a model $p(y|x,w)$ to predict the distribution of $y$ for a given new $x$, which requires parameters $w$. Which parameters to use? You can
- use a point estimate $w^*$ learned from the data $D$ (e.g. maximum likelihood estimate (MLE) or maximum a posteriori (MAP)). Predictions are given by $p(y|x,D) \approx p(y|x,w^*)$;
- or take a Bayesian approach by averaging over all possible values of $w$ weighted by $p(w|D)$, their probability given the observed data - a.k.a the posterior - learned from the data (e.g. by maximising the ELBO). Predictions are given by $p(y|x,D) = \int_w p(y|x,w)p(w|D)$.
Maximising the ELBO implicitly minimises $KL[q(w|D)\|p(w|D)]$ (see below). The ELBO is maximal when the KL term equals 0, meaning that the approximate posterior $q(w|D)$ equals the true posterior $p(w|D)$ (which is intractable to compute directly).
Longer answer:
Your definition of the ELBO isn't quite right.
You want to approximate the posterior $p(w|D)$ with $q(w)$ (whether this is denoted $q(w|D)$ or $q(w)$ doesn't particularly matter, but the former makes clearer it approximates the posterior, not a prior over $w$, $p(w)$, so I use $q(w|D)$ below).
To fit $q(w|D$ to $p(w|D)$, typically an intractable distribution, as you say, you maximise a lower bound (ELBO) on the log likelihood of the data, where $p(D) = \int_w p(D|w)p(w)$.
The ELBO os derived:
$\log p(D) = \int_w q(w|D)\log \tfrac{p(D,w)}{p(w|D)}
\quad= \int_w q(w|D)\log \tfrac{p(D|w)p(w)}{q(w|D)}\tfrac{q(w|D)}{p(w|D)}$
$= \int_w q(w|D)\log p(D|w) - \int_w q(w|D)\log \tfrac{q(w|D)}{p(w)} + \int_w q(w|D)\log \tfrac{q(w|D)}{p(w|D)}$
$\geq \int_w q(w|D)\log p(D|w) - \int_w q(w|D)\log \tfrac{q(w|D)}{p(w)} $
$= \mathbb{E}_{q(w|D)}[\log p(D|w)] - d_{KL}[q(w|D)\|p(w)]$
The difference to what you wrote is in the KL term. In particular the ELBO doesn't feature the true posterior $p(w|D)$ as you don't know it (if you did you wouldn't be trying to approximate it).
The dropped term leading to the inequality "gap" (leading to a lower bound), is the KL divergence between the approximate and true posterior, which is (implicitly) minimised as the ELBO is maximised. Subject to how well the various probabilities are modelled, the ELBO is fully maximised when the approximate posterior matches the true posterior
- (i.e. the KL term between the approximate and true posteriors goes to zero, hence the inequality "gap" disappears and you have an equality).
So maximising the ELBO is an elaborate indirect way of getting the approximate posterior $q(w|D)$ to fit the true posterior $p(w|D)$, which can then be used for predictions etc.
- It requires (i) an assumed prior over the parameters $p(w)$, (ii) the the data likelihood for a given parameter choice $p(D|w)$ to be computable, and (iii) a parametric family of distributions $q(w|D)$.
- Here, the likelihood $p(D|w) = \sum_{(x,y)\in D} p(y|x,w)$ is the same likelihood function as used for making predictions "at test time" (as you refer to at the outset).
- Given sufficient data, the posterior $q(w|D)$ may concentrate on a particular value $w^*$ (the maximum a posteriori estimate), allowing $w^*$ to be plugged directly into $p(y|x,w)$ for a good approximation to having to integrate over the whole posterior distribution.