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:
- calibrate the parameters using Bayesian inference, i.e., compute $p(\boldsymbol{\theta}|D)$
- 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:
- if the number of parameters $d$ is small enough, then numerical quadrature may be sufficient
- 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?