Implement inference for the latent field in INLA

[Michal-Novomestsky]

While inference on the hyperparameters $p(\theta \mid y)$ is relatively straightforward to obtain, efficiently estimating $p(x \mid y)$ is less trivial as the dimensionality of $x$ is in general very large.

Whilst Laplace approximations could either be applied on A. the joint $x$ or B. individually on each $x_i$, these approaches are too inaccurate and expensive respectively. Instead, a recommended middle-ground (used by R-INLA) is to use a bespoke method described by this paper. Note that it seems to natively make use of sparse tricks to work.

For a higher-level description, scroll down to "Bringing it all together: Integrated Nested Laplace Approximations" here.

In PyMC-Extras, I would recommend using get_laplace_approx for calculating the laplace approximations needed, however this may be insufficient for all of the bespoke steps described in the paper.