Inference: How is the Laplace approximation actually useful to us compared with MLE and MAP? I was reading a few different sources (including the "Machine Learning and Pattern Recognition" book by Bishop) about the Laplace integral approximation method for inference. However, I am still confused about the use of this method within inference.
Question: In what sense is the Laplace approximation better than MAP and MLE for inference calculations? Why would the enabling of the integration allow for better results?
Attempt: My understanding is as follows based on the reading:
In working probabilistically with parameters, we start by writing the probability of everything, also called the generative model. Then we get:
$$ p(f_{\star} | \mathcal{D}) = \frac{p(f_{\star}, \mathcal{D})}{p(\mathcal{D})} = \frac{\int p(f_{\star}, \mathcal{D}, \theta) d\theta}{p(\mathcal{D})} = \frac{\int p(f_{\star}| \mathcal{D}, \theta) p(\mathcal{D} | \theta) p(\theta) d\theta}{p(\mathcal{D})} $$
(However, I am not quite sure what $ f_{\star} $ is).
Then, the literature goes on to say that the core challenge of probabilistic inference is integration, marginalizing over the unknown, e.g
$$ p(f_{\star} | \mathcal{D}) = \int p(f_{\star} | \mathcal{D}, \theta) p(\theta | \mathcal{D}) d\theta $$
and
$$ p(\mathcal{D}) = \int p(\mathcal{D} | \theta) p(\theta) d \theta $$
Marginalizing can be interpreted as averaging over possible $\theta$ values, weighted by posterior probabilities (is this referring to the first or second integral?). Remember that maximum likelihood aimed to maximize $p(\mathcal{D} | \theta)$ and MAP aimed to maximize $p(\theta | \mathcal{D})$. These approximated the aforementioned expressions by delta functions to resolve tractable integrals. Then the reading says: "Better alternatives come up with more accurate models of those pdfs that nonetheless allow integration to be performed."
[EDIT #1]: I understand the mathematics of the method (i.e. using Taylor approximation to find mode of distribution we want to estimate; fitting a gaussian to approximate the function which then allows us to use standard results to compute the integral), but I just don't understand what the point of that is (basically a 'so what?' question)

Maximum likelihood is about maximizing the likelihood function to find the point estimate of the parameters
$$
\hat \theta_\text{ML} = \operatorname{arg\,max}_\theta \; p(X | \theta)
$$
In the case of Bayesian inference, we are looking for posterior distribution of the parameters, i.e. we consider also the prior $p(\theta)$
$$
p(\theta | X) = \frac{ p(X | \theta) \, p(\theta) }{ \int \, p(X | \theta) \, p(\theta) \,d\theta }
$$
the problematic part is usually calculating the integral in the denominator. Hopefully, we can approximate the posterior distribution by using Markov Chain Monte Carlo to sample from it and treat the empirical distribution of the samples as an approximation of the posterior. There are also other approximations possible. Maybe you don't need to know the full posterior distribution and knowing the mode of it would be enough? In such a case, you can use maximum a posteriori estimation
$$
\hat \theta_\text{MAP} = \operatorname{arg\,max}_\theta \; p(X | \theta) \, p(\theta)
$$
(We don't need the normalizing constant for optimization.)
Another approximation is Laplace approximation, where we approximate the posterior distribution with Gaussian distribution. In such a case, we get something in-between using MAP and full Bayesian inference, because we are left with distribution, but not an exact one, while it is just a little bit more complicated to do than using MAP alone.
