I have a nonlinear function of four variables $f(x_1,x_2,x_3,x_4)$. This function outputs three things: $[y_1, y_2, y_3]$.
(In reality, the function $f$ is a simulation which solves the equations of galactic evolution. The input parameters $x_1$ and $x_2$ are the mass and age of a galaxy, and $y_1$ and $y_2$ are its radius and luminosity at that particular time. The other inputs and outputs are more complicated to describe.)
This function is very expensive to compute, taking ample time on a computing cluster. I have evaluated this function with many different input values and stored a table of (inputs, outputs) pairs. The inputs were selected randomly from a pre-selected uniform interval.
I have measurements of an object (a real observed galaxy) which I would like to interpret in terms of my function $f$. In particular I have measured values $\tilde y_1, \tilde y_2, \tilde y_3$ and their associated (estimated) uncertainties $\sigma_1, \sigma_2, \sigma_3$.
(In other words, someone has measured the luminosity and some uncertainty for a galaxy, as well as some other measurements, and I would like to make a statement about the possible age and mass of the galaxy, assuming my simulations are correct.)
I want to know which values of $x_1, x_2, x_3, x_4$ and their uncertainties which are compatible with my measured values. I want to compute the credible region.
I can define a fitting function like $\sum_i^3 (y_i - \tilde y_i)^2/\sigma_i^2$, and combined with my priors on the $x$s (for example I know from previous experiments that $x_3$ is normally distributed with mean 5 and standard deviation 0.1), I can go through my table and find the best matches, and then look at the associated $x$ values.
But how do I estimate the credible region over the $x$s? Presumably, this region will depend on the resolution of the grid, and how well any of the other rows of the table compare to the 'best' one, etc. Is there a standard procedure?
I know how to solve this problem using MCMC. I can define my likelihood function as before and also encode my priors, and then churn on $f(\mathbf x)$. But now I have a table of values, and I do not trust interpolation in this table. So how can this calculation be made from a pre-computed grid?