Performing Bayesian prediction in practice when an explicit expression for the likelihood is not available

Question

I have a numerical code which, given an input vector $\mathbf{x}$ and parameter vector $\boldsymbol{\theta}$, gives me an output $l=f(\mathbf{x},\boldsymbol{\theta})$. I assume the statistical model

$y=f(\mathbf{x},\boldsymbol{\theta})+\epsilon$

where $\epsilon\sim\mathcal{N}(0,\sigma)$. I have a random sample $D=\{\mathbf{x}_i,y_i\}$ and I would like to:

1. calibrate the parameters using Bayesian inference, i.e., compute $p(\boldsymbol{\theta}|D)$
2. use $p(\boldsymbol{\theta}|D)$ to perform output predictions with associated credible intervals.

I know how to proceed in theory, but I'm not sure about the fiddly bits. So, as always

$p(\boldsymbol{\theta}|D)=\frac{p(D|\boldsymbol{\theta})p(\boldsymbol{\theta})}{\int p(D|\boldsymbol{\theta})p(\boldsymbol{\theta})d\boldsymbol{\theta}}$

The likelihood should be:

$p(D|\boldsymbol{\theta})=\frac{1}{\sqrt{2\pi}^N\sigma^N}\prod_{i=1}^N \exp{\left(-\frac{(y_i-{f_i(\mathbf{x}_i,\boldsymbol{\theta}))^2}}{2\sigma^2}\right)}$

Since the code is a black box, I don't have an explicit expression for $p(D|\boldsymbol{\theta})$, but I can compute this expression for any given $\boldsymbol{\theta}$. I know $p(\boldsymbol{\theta})$ (I choose it). The hard part is computing

$\int p(D|\boldsymbol{\theta})p(\boldsymbol{\theta})d\boldsymbol{\theta}$

I can imagine two approaches:

1. if the number of parameters $d$ is small enough, then numerical quadrature may be sufficient
2. otherwise, I think I should use an MCMC code. Suppose I have one available: the algorithm asks me in input the unnormalized target distribution (see for example here), which would be (correct me if I'm wrong) $p(D|\boldsymbol{\theta})p(\boldsymbol{\theta})$. The problem is that I don't have an explicit expression to pass to the algorithm: I can only evaluate it for given $\boldsymbol{\theta}$. So is the only solution to write my own MCMC? If so, I'll ask another question for details on how to do this.

Anyway, suppose I managed to compute $\int p(D|\boldsymbol{\theta})p(\boldsymbol{\theta})d\boldsymbol{\theta}$ one way or another. Now I have an "expression" for $p(\boldsymbol{\theta}|D)$, meaning that I can compute $p(\boldsymbol{\theta}|D)$ for any given $\boldsymbol{\theta}$. If I could sample from this distribution, then I think I could perform prediction this way: I choose a new input vector $\mathbf{x}^*$ where I want my prediction. I draw a sample of size $m$ $\{\boldsymbol{\theta}^{(1)},\dots,\boldsymbol{\theta}^{(m)}\}$ from $p(\boldsymbol{\theta}|D)$. I then compute the size $m$ sample

$\{y^{(1)}=f(\mathbf{x}^*,\boldsymbol{\theta}^{(1)}), \dots,y^{(m)}=f(\mathbf{x}^*,\boldsymbol{\theta}^{(m)})\}$

This would be my prediction at $\mathbf{x}^*$. Right? The problem is how to sample from $p(\boldsymbol{\theta}|D)$. What approaches are available to do this?

PS if $\sigma$ is unknown, I think I can still apply the same approach by just including it in $\boldsymbol{\theta}$. Correct?

Answer 1

The problem you have described, integrating $\int p(D|\boldsymbol{\theta})p(\boldsymbol{\theta})d\boldsymbol{\theta}$, is the driver of the notorious computational difficulties of Bayesian analysis. MCMC seeks to avoid evaluating this integral directly by sampling from the joint posterior $p(\boldsymbol{\theta}|D)$. As jaradniemi says, MCMC is not the only solution, but it is far and away the most common solution among Bayesian practitioners these days. A few other approaches in no particular order: importance sampling, quadrature, grid approximation, rejection sampling.

There are also various approximations for the posterior distribution, like fitting a multivariate normal (sometimes you'll see this kind of approximation justified through the so called "Bayesian CLT") or using variational bayes. The approximations have the advantage of speed, but it's often difficult to assess their accuracy without doing the full inference through some other method.

Assuming whatever method you use yields a sample $\{\boldsymbol{\theta}^{(1)},\dots,\boldsymbol{\theta}^{(m)}\}, \ i=1,\ldots,m$ from the joint posterior (as you would get from MCMC, for example), to handle a new input $\boldsymbol{x}^*$, you would draw $\{y^{(1)}, \ldots, y^{(m)}\}$ from $y^{(i)} \sim \operatorname{\mathcal{N}}\left( f(\boldsymbol{x}^*, \boldsymbol{\theta^{(i)}}),\sigma \right)$. This would give you the posterior preditive distribution for $y$ corresponding with the input vector $\boldsymbol{x^*}$. If desired, you could take whatever summary statistic you wanted from this distribution -- the mean, median, another quantile, etc.

Edit showing an example of grid approximation

Grid approximation is a simple idea where you evaluate an un-normalized distribution across a grid of points, and then gather samples from the grid points. For it to work well, you must have a fine enough grid. Here's 1-dimensional example from a normal distribution using R. As you can see, it's easy to extend to multiple dimensions, but the computation gets expensive quickly as the number of dimensions grows. For example, if you have 3-dimensions, even if your grid is only $1000$ points in each dimension, you need $1000^3$ grid points.

mu.truth <- 0                                                               
sd.truth <- 1.5                                                             

x.grid <- seq(-10, 10, length.out=10000)                            
some.constant <- 8585 # Totally arbitrary, just to illustrate we don't need a normalized density to use grid approximation
q.dens <- dnorm(x.grid, mean=mu.truth, sd=sd.truth)/some.constant       
normalized.dens <- q.dens/sum(q.dens) # R's sample function does this automatically, but just to be explicit...
y.sample <- sample(x.grid, 1000, replace=TRUE, prob=normalized.dens) # Gives us a sample from N(mu.truth, sd.truth^2)

# Visualize how well the approximation did                                  
hist(y.sample, freq=FALSE)                                                  
lines(x.grid, dnorm(x.grid, mean=mu.truth, sd=sd.truth), col="blue") 

As you can see, this results in a relatively representative sample from the un-normalized pdf.