Find shape and scale of InvGamma distribution

fitInvGamma.py

# how to do the same as Michael Betancourt in
# https://betanalpha.github.io/assets/case_studies/gp_part3/part3.html

import numpy as np
from scipy.stats import invgamma
from scipy.optimize import minimize

f = lambda x, lims: \
    (np.array([invgamma(a=x[0], scale=x[1]).cdf(lims[0]) - 0.01, 
               invgamma(a=x[0], scale=x[1]).sf(lims[1]) - 0.01])**2
    ).sum()


result = minimize(f, x0=[8, 30], args=[2, 10],
                  bounds=[(0, None),(0, None)], tol=1e-10)


a, b = result.x