How do I extract the maximum likelihood point estimates from a Turing output chain

Specific Domains Probabilistic Programming
1

I have a Turing output computed something like

model = make_model(data, ps)
chain = sample(model, NUTS(), MCMCSerial(), n_steps, n_chains; progress = false);
plot(chain)

image
image1009Ă—754 177 KB

How do I extract the maximum likelihood point estimate values, i.e. the peaks of the distributions? Is there some purposefully made function for doing this?

1 Like
2

Hi @Torkel !

There are two functions in Turing that will let you get the maximum likelihood or the maximum a posteriori.

# Generate a MLE estimate.
mle_estimate = maximum_likelihood(model)

# Generate a MAP estimate.
map_estimate = maximum_a_posteriori(model)

Notice that it is possible to use some arguments (the minimizer employed, the initial point for the minimization) to improve the results (but for a 3 dimensional problem, this should be straighforward :slight_smile: ).

1 Like
3

Thanks!

In this case this involves running the full MCMC chain though?I already have a (four) chain(s), so all the “heavy lifting” is already done. Do you know if there is a way to extract these values from the chain that have already been computed?

4

No, the mode estimation is (at least in the case I use it!) orders of magnitude faster:)
Also: notice that if you take the point of the chain with the highest value of the joint posterior, this will differ from the MLE in general (but this might depend on the prior you choose) and this is being a noisy estimate of the MAP as well. How many parameters in your model? If it’s just three, the mode estimation functionality should be pretty fast.

1 Like
5

Only 3, will try this, thanks! :slight_smile:

6

Is there a way to actually compute this directly from the chain? If not, I will mark the question as solved.

7

It could be possible to access informations related to che comptued chains (through chn.value.data), but yeah, for precise result I recomment to use the mode estimation utils (and it will be also faster than running chains).

8

If you have the chain already, you can easily find the sample with the highest log-probability. I’d do it this way:

df = DataFrame(chain)
df[argmax(df.lp), :]

As @marcobonici said, this won’t be the exact mode, but should be pretty close if you’ve collected enough samples. If your model is small and runs fast, no reason not to just use maximum_a_posteriori and get it exactly, even if you’ve already done MCMC.

3 Likes
9

Eventually, you can get this “slightly” suboptimal maximum a posteriori and use as an initial point for your minimization.

1 Like
10

Thanks both of you, very helpful!