# Computing the credible region from a pre-computed list of function evaluations

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?

• Build an emulator for f, then use MCMC for the unknown x values. But it concerns me that you say "I do not trust interpolation" which may indicate that f is not at all smooth which, in turn, means that your uncertainty in x will be large. – jaradniemi Apr 4 '19 at 17:04
• @jaradniemi The problem I am facing with interpolation is that the interpolation error is larger than any of my $\sigma$s (because the precision on the data are actually unprecedented). I can see this for example if I leave out an entry in the table and then try to predict its values using the interpolating function. – rhombidodecahedron Apr 4 '19 at 19:11
• @jaradniemi Another issue is that one of the properties, $y_3$, is actually undefined for some inputs of $\mathbf x$. I am not sure how to handle this case when interpolating. (Think something like: length of the bar of the galaxy. But some kinds of galaxies don't have a bar. I could omit the bar-less galaxies, but then it's a "hole" in parameter space.) – rhombidodecahedron Apr 4 '19 at 19:13