# Question about port of R code from the library “rethinking” to PyMC3

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
library(rethinking)

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


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.comJul 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.

• That std_q line is not in the R code - at least not in the part you quote. – Glen_b Jul 26 '18 at 1:29
• 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. – Glen_b Jul 27 '18 at 6:48