A very generous human named Osvaldo Martin did us the favor of porting all the R sample code in Richard McElreath's superb book Statistical Rethinking to PyMC3. I'm hugely grateful, but I've already encountered an example where the port depends on some knowledge about what the R code is doing "under the hood" and I would like to know what algorithm is being implemented.

Here's the R code (which uses McElreath's library "rethinking"):

# 2.6 - MAP

globe.qa <- map(
    w ~ dbinom(size = 9, prob = p),
    p ~ dunif(min = 0, max = 1)

which is ported to the following PyMC3 code:

data = np.repeat((0, 1), (3, 6))
with pm.Model() as normal_approximation:
    p = pm.Uniform('p', 0, 1)
    w = pm.Binomial('w', n=len(data), p=p, observed=data.sum())
    mean_q = pm.find_MAP()
    std_q = ((1/pm.find_hessian(mean_q, vars=[p]))**0.5)[0]
mean_q['p'], std_q

I am new to Bayesian statistics, but I know that Hessians are used in quadratic approximations, and I assume something of that sort is in play here, but where does this precise formula come from?

std_q = ((1/pm.find_hessian(mean_q, vars=[p]))**0.5)[0]

migrated from stats.meta.stackexchange.com Jul 25 '18 at 16:44

This question came from our discussion, support, and feature requests site for people interested in statistics, machine learning, data analysis, data mining, and data visualization.

  • 1
    $\begingroup$ That std_q line is not in the R code - at least not in the part you quote. $\endgroup$ – Glen_b Jul 26 '18 at 1:29
  • $\begingroup$ Even if it wasn't clear from the question, it should be reasonably clear from the answer that this is "a question requiring statistical expertise to answer", and thereby clearly on topic. I will reopen on that basis; please feel free to take issue with it on meta. $\endgroup$ – Glen_b Jul 27 '18 at 6:48

The Hessian matrix is widely used in statistics to obtain an asymptotic approximation to the covariance matrix of the parameter estimates in large sample problems (specifically from its inverse). I presume this particular calculation is in a Bayesian context in which case its not exactly that - it would be an approximation to the covariance matrix of the joint conditional posterior distribution of parameters.

In turn that covariance matrix is mainly used to pull out estimated standard errors (standard deviation of the marginal posterior distribution of the parameter conditional on the data in a Bayesian context), by taking the square roots of diagonal elements.

  • $\begingroup$ Thanks very much for taking time to answer, Glen. The information you provided led me to a tutorial on Fisher Matrices which together with your summary has helped me get an idea of what's going on in that code. Here's a link for anyone else who's interested: youtube.com/watch?v=m62I5_ow3O8 $\endgroup$ – John Strong Jul 26 '18 at 12:14
  • 1
    $\begingroup$ By the way, Glen, you correctly pointed out that the line of code assigning a value to std_q is not in the R code, but that's what I meant when I said the map() method from the rethinking library was generating an STDEV "under the hood". I assume the person who ported the map() method to PyMC3 must have viewed the R code in the map() function to figure out what it was doing, and he discovered something like an inverse Hessian. It is not clear from the R code snippet, but map() returns two values: a mean and a STDEV, just like the PyMC3 code. $\endgroup$ – John Strong Jul 27 '18 at 2:58

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.