More flexible fitting function, allow likelihood, remove uncertainties dependency

[aminnj]

Based on the discussion in #146, I added some features to plot_pull. The theme is making fitting more streamlined for exploration.

• Pull curve_fit initial guess (p0) from default arguments, if they exist
• Allow a string as an alternative to a lambda function (plot_pull("a+b*x"))
• Cosmetic change to the band, and embedding fit result into the legend
• Likelihood fit (plot_pull(..., likelihood=True)) (chi2 by default, as before)
• Remove uncertainties.numpy dependency and construct band by resampling covariance matrix
• Introduce iminuit (gets the covariance matrix right, unlike scipy.optimize most of the time), but the initial guesses for iminuit are seeded from scipy

Setup

import numpy as np
from hist import Hist

np.random.seed(42)
hh = Hist.new.Reg(50, -5, 5).Double().fill(np.random.normal(0,1,int(1e5)))

Before (including a bug-fix for the variances from the above issue):

from uncertainties import unumpy as unp
def func(x, constant, mean, sigma):
    exp = unp.exp if constant.dtype == np.dtype("O") else np.exp
    return constant * exp(-((x - mean) ** 2.0) / (2 * sigma ** 2))

hh.plot_pull(func)

After:

# as before, but no need for `uncertainties.numpy` as the error band comes
# from resampling the covariance matrix
def func(x, constant, mean, sigma):
    return constant * np.exp(-((x - mean) ** 2.0) / (2 * sigma ** 2))
hh.plot_pull(func)

# `curve_fit` `p0` extracted from defaults, if any
def func(x, constant=80, mean=0., sigma=1.):
    return constant * np.exp(-((x - mean) ** 2.0) / (2 * sigma ** 2))
hh.plot_pull(func)

# strings are allowed to allow for more compactness than a lambda
# x is assumed to be the main variable
hh.plot_pull("constant*np.exp(-(x-mean)**2. / (2*sigma**2))")

# gaussian is a common/special function, so this also works
# reasonable guesses are made for constant/mean/sigma
hh.plot_pull("gaus")

# chi2 puts `a` around 5, but likelihood puts `a` around 1e3/50 = 20
hh.plot_pull("a+b*x", likelihood=True)

[henryiii]

Can you verify my changes are valid? Could you add a test with pytest-mpl, perhaps?

[henryiii]

Notes on my changes.

[henryiii]

This is the default for Callable - best to not leave empty Generics.

[henryiii]

ArrayLike can be a simple number, so xdata[stuff] was not valid. Either use np.asarray, or just require arrays.

[henryiii]

Makes the typing easier, because NumPy's typing is a bit spotty. It doesn't know that *array is valid, etc. This is small so was a simple fix.

[henryiii]

I'm using inspect.signature, as it's a more public interface; __defaults__ (and maybe __code__) are untyped, AFAICT (or maybe just not always present on Callables). This also supports partial defaults, while the other was all or nothing.

[henryiii]

This could be a ModuleNotFoundError.

[henryiii]

Never check the type exactly, always use isinstance. It supports subclassing and MyPy type narrowing.

[henryiii]

Also never make a list for containership testing, use a set, it's faster. x in {a, b}.

[henryiii]

Not always fond of NumPy's choices here. Not sure what a number[Any] is. So just forcing it to a Python float.

[henryiii]

Same higher level usage of inspect.signature over .__code__.co_varname.

[henryiii]

@all-contributors please add @aminnj for code

[allcontributors]

@henryiii

I've put up a pull request to add @aminnj! 🎉

[henryiii]

Also, would you like to fill in the changelog so we can clear that needs-changelog badge? :)

(I can do it if you prefer)

[aminnj]

I appreciate the comments. I haven't yet grasped all the subtleties of the type system...

Is the suggestion to use pytest-mpl to make a snapshot of an example plot and use it as a reference image in unit testing? What if the style is changed upstream?

[henryiii]

This also need docs, but don't worry about that - I'll be reworking the docs in the near future to make them easier to edit. When I do, I'll add to this part.

There are two ways to test. The easier way is to use pytest-mpl to make a snapshot of an example plot and use it as a reference image in unit testing. If the style is changed upstream, we will immediately know and will have to change our test image.

The better way to test is to mock mpl and then verify the sequence of commands to the mpl functions do not change. That also gives much better error messages if something gets changed, rather than just an image diff. But that's much harder to set up. You can see what I did for mplhep here: https://github.com/scikit-hep/mplhep/blob/5effa0b857d7eea93da0299cbbcdc777bd3abaaf/tests/test_mock.py - but as you see, I didn't end up adding it for everything because it's a bit of work for each one.

[aminnj]

I added in pytest-mpl and seems CI is finally happy with it. Lowered the DPI for the test case made with pytest tests/test_plot.py --mpl-generate-path=baseline. By default, it has a very matplotlib 1.x feel.

[henryiii]

Looks good to me, can you add a changelog entry, then we can merge? If not, let me know, and I'll add one; I usually lightly push contributors to add changelog entries but add them if not; original contributors are usually the best at advertising what they have done. :)

[aminnj]

I added in some test cases that I had forgotten earlier (the h.plot_pull("gaus") alias and likelihood=True). The changelog entry I added is just a compact summary of the points in the original post of this thread, but let me know if I should include more details.

[henryiii]

Looks great, thanks!

[eduardo-rodrigues]

Hi @henryiii, @aminnj, I'm following your nice developments. I just wanted to make a comment now that I see h.plot_pull("gaus"): I know that ROOT does this, amputating Gauss's name to save a little letter, but I hate that. TBH I see no gain and consider far more important to save the great mathemacian's name. Could we please please do "gauss"?

[aminnj]

I personally have no strong attachment to "gaus", so I can offer a +1 to "gauss"

[nsmith-]

I was taking a look at https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.html today and noticed it has a fit function that does an unbinned max-likelihood fit. I wonder if there's a possibility we could add a binned_fit function to scipy? It could take a UHI plottable or numpy array as input.