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

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?

• It's not clear, to me, what you can do with $f$ so you may want to give more background on $f$. You state "I can ... evaluate [the unnormalized target distribution] for a given $\theta$." This is what is necessary to perform most MCMC so yes, MCMC seems like the way to go although it is never the "only" solution. – jaradniemi Jan 24 '17 at 12:31
• @jaradniemi thanks for the comment. Consider you have a large computational code, for example a FEM code to perform structural analysis, a CFD code for the numerical solution of fluid flows, a meteorological code to predict the temperature tomorrow at your place, etc.. You can run the code, meaning that you can give it an input vector $\mathbf{x}$ and a calibration parameter vector $\boldsymbol{\theta}$. The code will return you an output. Since I'm thinking of a deterministic code, each time I run the code with the same inputs and parameters I will get the same output. This means that I [1\2] – DeltaIV Jan 24 '17 at 12:41
• [2/2] can consider the code as a function $f(\mathbf{x},\boldsymbol{\theta})$. I don't have the explicit expression of the function, but for each couple $(\mathbf{x},\boldsymbol{\theta})$ I can run the code and get the output $l=f(\mathbf{x},\boldsymbol{\theta})$. This is all I know about $f$. – DeltaIV Jan 24 '17 at 12:44
• @DeltaIV. You do not need an explicit expression for the function and this is not so uncommon (just think about truncated normal distribution). I think your question is more about finding a MCMC library allowing you to provide a sophisticated likelihood function. In this case, it would be helpful to give details about your programming environment. – peuhp Jan 24 '17 at 13:00
• @peuhp R and MATLAB. I can use Python if the need arises, but I'd really prefer to use the other two. So do you think that the rest of my approach is sound? – DeltaIV Jan 24 '17 at 13:14

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.

• +1. I don't know about grid approximation (a link would be nice). I know about quadrature, but quadrature doesn't give me a sample $\{\boldsymbol{\theta}^{(1)},\dots,\boldsymbol{\theta}^{(m)}\}, \ i=1,\ldots,m$. It only allows me to compute $p (\boldsymbol{\theta}|D)$ for any $\boldsymbol{\theta}$. How can I use this to sample from $p(\boldsymbol{\theta}|D)$? And if I use MCMC instead, which gives me such a sample, can I sample with replacement from the MCMC sample to get new $y$ samples? I think so because this would basically correspond to bootstrap. – DeltaIV Jan 26 '17 at 7:11
• Yes, quadrature won't give you a sample. To get a sample you'll need to use something like rejection sampling or importance sampling. I wasn't sure about a good link for grid approximation, so I added a little example. It's covered in many texts though -- I know for sure grid approximation is covered in Gelman's book and in Kruschke's book. I don't think it really makes sense to sample with replacement from an MCMC sample. You get full posterior inference from MCMC, so I can't think of any reason to do a bootstrap. Maybe there's some part of your problem I'm not understanding though. – AtALoss Jan 27 '17 at 7:42