use MCMC for utils.get_walkers positions

[kecnry]

Is your feature request related to a problem? Please describe.

Using PRISM to probe a large parameter space and then calling utils.get_walkers can be very expensive as the whole parameter space needs to be re-evaluated for a small percentage of plausible regions.

Describe the solution you would like

Allow for an alternate option to start at pipe.impl_sam and explore nearby parameter space until hitting implausible regions.

[1313e]

@kecnry Yeah, finding a large number of plausible samples in a small plausible space can be very expensive. For analyzing the emulator, it is required to use Latin-Hypercube sampling in order to avoid biasing. But, for obtaining plausible MCMC starting positions, an alternative way would be possible.

I will have a look at that and see if I can implement such a thing in a nice way.

[1313e]

@kecnry So, some updates: I got a somewhat working implementation right now, but it is not fast enough when dealing with very small parameter spaces. As I have a few very important deadlines coming up, I am not sure how soon I can do the proper implementation. It might take until the end of next week or even longer.

Edit: Never mind, I think I have found an implementation that is pretty quick.

[1313e]

@kecnry I have now implemented the new functionality into the 'utils.get_walkers' function. Can you check if it performs as you expect?

[kecnry]

It seems to work! The example we discussed that would have taken days now takes 20-30 seconds to get the walkers, so that's definitely promising.

Just to be clear: if I don't pass init_walkers, then this starts at pipe.impl_sam and otherwise it'll re-search the whole grid with roughly-evenly-sampled init_walkers (so essentially just doing Metropolis-Hastings over the whole parameter space)? I read the docstring, but couldn't really figure out what the expected behavior should be in this case.

In the cases where 0 plausible samples are returned by analyze, would you suggest using this with init_walkers or re-running with a more-conservative impl_cut or finer sampling?

[1313e]

Awesome. And no, both input arguments do something differently. If you give a value to init_walkers, it will create an LHD of that size and evaluate it in the emulator, like before. Not giving it will use pipe.impl_sam. If you also give a value for req_n_walkers, then these will be used as the starting points for all MCMC chains in the MH sampling.

So, the way that the init_walkers argument works has not changed. What is different now is that when you provide a value for req_n_walkers, then instead of returning all plausible samples in init_walkers, it uses them to find the required number of plausible starting positions.

And in the case that pipe.analyze returns no plausible samples, I would first check what the residual variance is (pipe.emulator.rsdl_var). If they are very low, you need finer sampling. If not, you need more conservative cut-offs. You always have to make sure that pipe.analyze can return plausible samples, as it is the only method in the pipeline that is allowed to modify the implausibility parameters. If it cannot return plausible samples, then get_walkers probably will not either, regardless of how many samples you throw in.

[1313e]

You do have to be careful though. If the plausible samples in init_walkers do not cover all plausible regions because they were too sparsely sampled, the MH sampling algorithm will not find them either. An example of this can be seen in the following:

The top plot is the projection figure of an emulator that I was testing the algorithm with (GaussianLink). I then used 7 plausible samples to find 10k plausible starting positions, which are shown in the bottom plot. As you can see there, it captures the main plausible region perfectly fine (even the density), but there was still a very small plausible region at the top left (around [1.5, 9.0]). Although it is possible that sampling more densely for the projections would have removed that region, it is currently a plausible region that was not covered with the MH sampling.