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?

  • $\begingroup$ 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. $\endgroup$ – jaradniemi Apr 4 '19 at 17:04
  • $\begingroup$ @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. $\endgroup$ – rhombidodecahedron Apr 4 '19 at 19:11
  • $\begingroup$ @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.) $\endgroup$ – rhombidodecahedron Apr 4 '19 at 19:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.