search

scikit-hep / iminuit

Jupyter-friendly Python interface for C++ MINUIT2
https://scikit-hep.org/iminuit
Other
283 stars 76 forks source link

Pass `ncall` and `iterate` in methods that build profile likelihoods #1008

Closed ikrommyd closed 1 month ago

ikrommyd commented 1 month ago

Hi @HDembinski,

Would you be open to allow migrad optional arguments as arguments to method that return profile likelihoods like mnprofile. I'm getting the warning that migrad failed to converge quite often in such cases and I thought this might be helpful. I could make a PR, just asking if you'd be open to it.

ikrommyd commented 1 month ago

In general, I think a little bit more customizability of the minimization during profiling doesn't hurt.

In my use case, I had to hack iminuit like this to get better convergence during scans:

diff --git a/src/iminuit/minuit.py b/src/iminuit/minuit.py
index 4e40958..889ad89 100644
--- a/src/iminuit/minuit.py
+++ b/src/iminuit/minuit.py
@@ -1623,6 +1623,12 @@ class Minuit:
             state.set_value(ipar, v)
             migrad = MnMigrad(self._fcn, state, self.strategy)
             fm = migrad(0, self._tolerance)
+            if not fm.is_valid:
+                warnings.warn(
+                    "Running simplex and then migrad again"
+                )
+                fm = MnSimplex(self._fcn, fm.state, self.strategy)(0, self._tolerance)
+                fm = MnMigrad(self._fcn, fm.state, self.strategy)(0, self._tolerance)
             if not fm.is_valid:
                 warnings.warn(
                     f"MIGRAD fails to converge for {pname}={v}", mutil.IMinuitWarning
HDembinski commented 1 month ago

It is a good idea, thanks. Will do.

ikrommyd commented 1 month ago

I can write up a PR, I don't know how much customizability you'd want though. I think ncall and iterate as arguments as well an option to use migrad or simplex as the minimizer is enough. If I user wants more custom stuff like a sequence of minimization algorithms etc, I think they can hack it themselves.

HDembinski commented 1 month ago

I rather do it myself, to keep the fitting consistent in various utilities. See the linked PR.

I want you to provide a minimal example in which scans fail in the current version without your hack. I can then check whether the update handles this.

ikrommyd commented 1 month ago

I don't know how minimal I can go when it comes to the model though. I've encountered this in a 1D line fit but the model is has different astrophysical signal and background components which are constructed from interpolating functions and the parameters are normalizations and spectral indices that change those components like component * energy^(spectral index). What I observed was that when doing the minos profile for one such normalization parameter, I would get a lot of warnings and points where the chi2 loss increased for no reason. After the hack, those disappeared. I tried some different minimization schemes but what makes the minos profile more robust is just simplex followed by migrad of course. Also I lowered the tolerance for the entire Minuit object and set the strategy to 2

ikrommyd commented 1 month ago

Now all that may also have to do with the fact that my model is not perfect in the current state. My parameters are bounded and during the fits I see a lot of parameters being at limit and having uncertainties as large as their allowed region (i.e. hitting an upper limit and having a 1 sigma minos error at the lower limit). I also performed a bayesian analysis to check this and sample the posterior and I see some 1D posteriors being flat or truncated by the boundaries. So I either need to rethink the modelling there or add some constraints to those parameters from the literature.

HDembinski commented 1 month ago

I understand it is often difficult to provide concrete examples in which iminuit fails.

If you cannot provide the example, then do you think you can install the iminuit version in the PR and see whether it works for you? There is a proper way and a shortcut, I would recommend the shortcut.

The shortcut is to install the latest release that you already have. Then find out where it is installed.

import iminuit
print(iminuit.__file__)

Note down the folder which contains the __init__.py. Checkout the PR and copy all the files in src/iminuit to this folder. Then try the updated version on your problem and report back whether it works as good as your hack, please.

ikrommyd commented 1 month ago

I got iminuit in editable mode and this is the code I'm running. I'll try later to provide a way more minimal example by pickling all the interpolating functions and providing the data and the model in a much simpler format.

import numpy as np
import numba as nb
import iminuit
from scipy.interpolate import interp1d
from matplotlib import pyplot as plt
from iminuit import Minuit
from iminuit.cost import LeastSquares

def convert_sim_file(path):
    df = np.loadtxt(path)
    E = df[:, 0]
    dNdE = df[:, 1]
    flux = E**2 * dNdE * 1e-4
    interp = interp1d(E, flux, kind="cubic")
    return interp(x), flux

def convert_sim_file(path):
    df = np.loadtxt(path)
    E = df[:, 0]
    dNdE = df[:, 1]
    flux = E**2 * dNdE * 1e-4
    interp = interp1d(E, flux, kind="cubic")
    return interp(x), flux

def read_sims(simtype, mass, channel, colossus_model, interp=True):
    fluxSFG = convert_sim_file(
        f"../FinalFilesForFitting_v2/StarForming_and_Starburst_Galaxies_With_Attenuation_{simtype}_v2.dat"
    )
    fluxRG = convert_sim_file(
        f"../FinalFilesForFitting_v2/RadioGalaxies_With_Attenuation_Tim_{simtype}_Model_v2.dat"
    )
    fluxFSRQ = convert_sim_file(
        f"../FinalFilesForFitting_v2/FSRQ_With_Attenuation_{simtype}_Model_IGRB_estimate_omega_0p9_v2.dat"
    )
    fluxBLLac = convert_sim_file(
        f"../FinalFilesForFitting_v2/BLLac_With_Attenuation_{simtype}_Model_IGRB_estimate_omega_0p9_mustar_{bllacversion}_v2.dat"
    )
    fluxUHECR = convert_sim_file(
        f"../FinalFilesForFitting_v2/UHECR_spectrum_protons.dat"
    )
    fluxMSP = convert_sim_file(f"../FinalFilesForFitting_v2/MSP_spectrum.dat")
    fluxDMEGRB = convert_sim_file(
        f"../DM_Intensity_Files/GammaRayIntensity_From_DM_annihilation_Fiducial_{channel}_{mass}GeV_sigmav_3e-26_{colossus_model}_HMF.dat"
    )
    fluxDMHALO = convert_sim_file(
        f"../DM_Intensity_Files/MilkyWayHalo_GammaRayIntensity_From_DM_annihilation_Fiducial_{channel}_{mass}GeV_sigmav_3e-26.dat"
    )
    if interp:
        ind = 0
    else:
        ind = 1
    return (
        fluxSFG[ind],
        fluxRG[ind],
        fluxFSRQ[ind],
        fluxBLLac[ind],
        fluxUHECR[ind],
        fluxMSP[ind],
        fluxDMEGRB[ind],
        fluxDMHALO[ind],
    )

def model(E, theta):
    NormSFG, DgammaSFG = theta[0:2]
    NormRG, DgammaRG = theta[2:4]
    NormFSRQ, DgammaFSRQ = theta[4:6]
    NormBLLac, DgammaBLLac = theta[6:8]
    NormUHECR, NormMSP = theta[8:10]
    NormDM = theta[10]

    sfg = NormSFG * fluxSFG * E ** (DgammaSFG)
    rg = NormRG * fluxRG * E ** (DgammaRG)
    fsrq = NormFSRQ * fluxFSRQ * E ** (DgammaFSRQ)
    bllac = NormBLLac * fluxBLLac * E ** (DgammaBLLac)
    uhecr = NormUHECR * fluxUHECR
    msp = NormMSP * fluxMSP
    dm = NormDM * (fluxDMEGRB + fluxDMHALO) * 1e4

    model = sfg + rg + fsrq + bllac + uhecr + msp + dm

    return model

df = np.loadtxt(f"../IGRB_Fermi_Spectrum_Model_A.dat")
simtype = "Fiducial"
bllacversion = "2p3"

x = df[:, 0]
y = df[:, 1]
upper = df[:, 2] - df[:, 1]
lower = df[:, 1] - df[:, 3]
yerr = np.sqrt(upper * lower)
yerr[-1] = upper[-1]

parrange = [
    [0.1, 1.5],
    [-0.3, 0.3],
    [0.2, 1.5],
    [-0.3, 0.2],
    [0.5, 1.5],
    [-0.04, 0.03],
    [0.2, 1.2],
    [-0.15, 0.15],
    [0.2, 5.0],
    [0.5, 1.5],
    [0, None],
]
theta0 = np.array([0.5, 0, 0.5, 0, 0.75, 0, 1, 0, 1, 1, 1])

mass = "10."
channel = "bbbar"
colossus_model = "Comparat17"
(
    fluxSFG,
    fluxRG,
    fluxFSRQ,
    fluxBLLac,
    fluxUHECR,
    fluxMSP,
    fluxDMEGRB,
    fluxDMHALO,
) = read_sims(simtype, mass, channel, colossus_model, interp=True)

jitted_model = nb.njit(model)
chi2 = LeastSquares(x, y, yerr, jitted_model)
m = Minuit(
    chi2,
    theta0,
    name=(
        "Norm SFG",
        "gamma SFG",
        "Norm RG",
        "gamma RG",
        "Norm FSRQ",
        "gamma FSRQ",
        "Norm BLLac",
        "gamma BLLac",
        "Norm UHECR",
        "Norm MSP",
        "Norm DM",
    ),
)
m.limits = parrange
m.strategy = 2
m.migrad()
m.hesse()
m.minos()
profile = m.mnprofile(
    "Norm DM", grid=np.logspace(-3 - np.log10(3), 3 - np.log10(3), 600)
)

With iminuit from the develop branch, I'm getting a lot of of warnings like:

[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.01012393923758175
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.011102385626231472
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.012459476070843946
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.0146425475078723
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.01887123318889198
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.027294012934681654
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.035176370947343175
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.03947612403608427
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.061185948976146015
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.07705816184989164
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.09267278210152602
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.10891032711910645
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.1371626942901648
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=1.020207303110914
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=77.95195155142515
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=79.77075212470125
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=83.53665398149796
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=93.74768417894836
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=159.34723118242457
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=166.86986344803555
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=170.76332598824624
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=174.74763207706894
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=178.82490130605981
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=182.9973027217824
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=187.267055979708
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=191.6364325250396
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=196.10775680108662
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=200.68340748583444
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=205.3658187573655
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=210.15748158880595
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=215.06094507348627
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=220.0788177810209
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=225.21376914502852
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=230.46853088323107
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=235.8458984506869
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=241.3487325269312
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=246.979960537815
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=252.7425782128522
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=258.639651178903
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=264.67431659104204
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=270.8497848014783
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=277.1693410674128
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=283.63634729875014
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=290.25424384658163
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=297.0265513334021
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=303.95687252602835
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=311.0488942522175
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=318.3063893620046
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=325.73321873480273
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): IMinuitWarning: MIGRAD fails to converge for Norm DM=333.3333333333332
  warnings.warn(

With the hack, I'm just getting the hack warning

[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1627](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1626): UserWarning: Running simplex and then migrad again
  warnings.warn(

but then the minima duing the scan are valid. With iminuit from the robuster_fitting branch, I'm still getting warnings, but a little bit less,

[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.07026707363549274
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.07358431580793595
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.09483505393512022
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.13403533770441942
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=0.14363801552448263
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=79.77075212470125
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=174.74763207706894
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=187.267055979708
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=230.46853088323107
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=235.8458984506869
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=241.3487325269312
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=246.979960537815
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=252.7425782128522
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=264.67431659104204
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=270.8497848014783
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=277.1693410674128
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=283.63634729875014
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=290.25424384658163
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=297.0265513334021
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=303.95687252602835
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=311.0488942522175
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=318.3063893620046
  warnings.warn(
[/Users/iason/fun/iminuit/src/iminuit/minuit.py:1651](http://localhost:8888/Users/iason/fun/iminuit/src/iminuit/minuit.py#line=1650): IMinuitWarning: MIGRAD fails to converge for Norm DM=325.73321873480273
  warnings.warn(
HDembinski commented 1 month ago

Thank you for testing this. It is strange that the new version differs from your hack, because they should be equivalent. The new code iterates MnSimplex and MnMigrad more than once, but only if the first iteration does not produce a valid fit. In your case it looks like the fit succeeds after one iteration.

It could be related to you using strategy 2. In these scans, I am reducing the strategy to 1 if you use 2 and to 0 if you use 1, because the higher strategies compute the Hesse matrix explicitly while iterating and that is very expensive.

However, in your case that maybe makes the difference, so I have to rethink this. Perhaps a good heuristic is to start with the faster strategy and fall back to a higher strategy automatically if the first iteration of Migrad fails.

ikrommyd commented 1 month ago

This is a version that loads the data and I've put the data and the interpolating functions in a pickle file in the attached zip

import pickle
import numpy as np
import numba as nb
import iminuit
from scipy.interpolate import interp1d
from matplotlib import pyplot as plt
from iminuit import Minuit
from iminuit.cost import LeastSquares

def load_data(filename):
    data = np.loadtxt(filename)
    return data[:, 0], data[:, 1], data[:, 2]

def load_interpolating_functions(filename):
    with open(filename, 'rb') as f:
        return pickle.load(f)

def model(E, theta):
    NormSFG, DgammaSFG = theta[0:2]
    NormRG, DgammaRG = theta[2:4]
    NormFSRQ, DgammaFSRQ = theta[4:6]
    NormBLLac, DgammaBLLac = theta[6:8]
    NormUHECR, NormMSP = theta[8:10]
    NormDM = theta[10]

    sfg = NormSFG * fluxSFG * E ** (DgammaSFG)
    rg = NormRG * fluxRG * E ** (DgammaRG)
    fsrq = NormFSRQ * fluxFSRQ * E ** (DgammaFSRQ)
    bllac = NormBLLac * fluxBLLac * E ** (DgammaBLLac)
    uhecr = NormUHECR * fluxUHECR
    msp = NormMSP * fluxMSP
    dm = NormDM * (fluxDMEGRB + fluxDMHALO) * 1e4

    model = sfg + rg + fsrq + bllac + uhecr + msp + dm

    return model

# Load data
x, y, yerr = load_data('fermi_spectrum_data.txt')

# Load interpolating functions
interpolating_functions = load_interpolating_functions('interpolating_functions.pkl')

# Assign the global variables
fluxSFG = interpolating_functions['fluxSFG'](x)
fluxRG = interpolating_functions['fluxRG'](x)
fluxFSRQ = interpolating_functions['fluxFSRQ'](x)
fluxBLLac = interpolating_functions['fluxBLLac'](x)
fluxUHECR = interpolating_functions['fluxUHECR'](x)
fluxMSP = interpolating_functions['fluxMSP'](x)
fluxDMEGRB = interpolating_functions['fluxDMEGRB'](x)
fluxDMHALO = interpolating_functions['fluxDMHALO'](x)

# Prepare model for fitting
parrange = [
    [0.1, 1.5],
    [-0.3, 0.3],
    [0.2, 1.5],
    [-0.3, 0.2],
    [0.5, 1.5],
    [-0.04, 0.03],
    [0.2, 1.2],
    [-0.15, 0.15],
    [0.2, 5.0],
    [0.5, 1.5],
    [0, None],
]
theta0 = np.array([0.5, 0, 0.5, 0, 0.75, 0, 1, 0, 1, 1, 1])

# JIT compile the model function
jitted_model = nb.njit(model)

# Fit the model
chi2 = LeastSquares(x, y, yerr, jitted_model)
m = Minuit(
    chi2,
    theta0,
    name=(
        "Norm SFG",
        "gamma SFG",
        "Norm RG",
        "gamma RG",
        "Norm FSRQ",
        "gamma FSRQ",
        "Norm BLLac",
        "gamma BLLac",
        "Norm UHECR",
        "Norm MSP",
        "Norm DM",
    ),
)
m.limits = parrange
m.strategy = 2
m.migrad()
m.hesse()
m.minos()

# Profile likelihood for "Norm DM"
profile = m.mnprofile(
    "Norm DM", grid=np.logspace(-3 - np.log10(3), 3 - np.log10(3), 600)
)
m

data_and_interpolating_functions.zip

HDembinski commented 1 month ago

To comment on your observation of the Bayesian posteriors: If you get flat posteriors that is indeed a sign that your fit has parameters which are not constrained by the data. That will make Minuit fail, and Minuit is not smart enough to detect this automatically and warn about the problem.

HDembinski commented 1 month ago

When I try your example (thanks), it fails with an error in the step m.minos() because the minimum is not valid.

ikrommyd commented 1 month ago

When I try your example (thanks), it fails with an error in the step m.minos() because the minimum is not valid.

Hmm, on my laptop I'm getting a valid minimum repeatedly (there's no RNG as well), but yes I've had this also for different fits though. I'm doing like thousands of such fits.

ikrommyd commented 1 month ago

To comment on your observation of the Bayesian posteriors: If you get flat posteriors that is indeed a sign that your fit has parameters which are not constrained by the data. That will make Minuit fail, and Minuit is not smart enough to detect this automatically and warn about the problem.

Yes I did this for debugging mainly, but even though I'm interesting in the upper limit on the dark matter normalization, I should probably constrain the astrophysical background more rather than having flat bounded ranges.

HDembinski commented 1 month ago

I have an Intel CPU, if you have a modern Mac you may have a M1 chip. This or other differences in used libraries can meddle with the results in pathologicial cases like this one.

ikrommyd commented 1 month ago

Yup, I'm on arm64 macos.

HDembinski commented 1 month ago

Anyway, I get the warnings, so I can try to debug it, thanks for now.

ikrommyd commented 1 month ago

The whole point of this issue was to basically just say that some more robustness and customizabilty would be nice. It's not to debug my exact and probably extreme example of course. I also don't know if I like the error when trying to run minos on an invalid minimum although I agree it makes sense from a statistics perspective. Don't know if Minuit2 or Minuit from ROOT gives a similar error or just runs minos anyways.

HDembinski commented 1 month ago

The error is thrown by iminuit, not ROOT. Calling Minos on an invalid minimum makes no sense. The algorithm assumes that it is in a valid minimum.

Since I don't want to use the hack you posted above in iminuit, I need a way to test whether my idea of handling this works as good or better. To check this I needed your example.

HDembinski commented 1 month ago

I did a few variations but cannot get rid of the fit failures with my code. Please submit your hack as a PR, so that I can test it on your example on my computer as well. It might be that I am getting the same failures. At that point I would give up trying to fix this.

The patch has the customizations that you want.

HDembinski commented 1 month ago

Looking at your concrete example again, the covariance matrix shows strong anti-correlations of Norm SFG, Norm RG, gamma RG, and as you say many parameters are at their limit. It looks like a pathological fit, with the model being underconstrained by the data.

ikrommyd commented 1 month ago

I did a few variations but cannot get rid of the fit failures with my code. Please submit your hack as a PR, so that I can test it on your example on my computer as well. It might be that I am getting the same failures. At that point I would give up trying to fix this.

The patch has the customizations that you want.

https://github.com/scikit-hep/iminuit/pull/1010 there you go

ikrommyd commented 1 month ago

The error is thrown by iminuit, not ROOT. Calling Minos on an invalid minimum makes no sense. The algorithm assumes that it is in a valid minimum.

Since I don't want to use the hack you posted above in iminuit, I need a way to test whether my idea of handling this works as good or better. To check this I needed your example.

Yup, I know the error comes from python, I just didn't know if the C++ Minuit code cared at all about the validity of the minimum.

And yes, of course the hack makes no sense to have in the central iminuit code.

HDembinski commented 1 month ago

The C++ code generally does not care a lot about internal consistency.

ikrommyd commented 1 month ago

Looking at your concrete example again, the covariance matrix shows strong anti-correlations of Norm SFG, Norm RG, gamma RG, and as you say many parameters are at their limit. It looks like a pathological fit, with the model being underconstrained by the data.

If you also look at this example fit, it looks like a lot of the background components can be used interchangably to make the fit work. You can lower the normalization of one and increase the normalization of another and not a lot will change in the combined line. I assume that's why I'm getting these "at limit" parameters and flat or truncated posteriors in the bayesian case.

image
HDembinski commented 1 month ago

You got it, that is the issue.

HDembinski commented 1 month ago

I wrote a new tutorial based on our discussion, which lists the most common reasons why fits fail or are unstable. There are some formatting mistakes in there that I will fix later, but I think it is a good read.

https://scikit-hep.org/iminuit/notebooks/unstable_fit.html

The latest release 2.27 has the changes that we discussed here

ikrommyd commented 1 month ago

Thank you Hans for the PR and the tutorial!!. Should be helpful to people struggling with such fits.

ikrommyd commented 1 month ago

You got it, that is the issue.

And you can also see why only the Ultra high energy cosmic-ray (UHECR) component is constrained well and also has a good-looking posterior. The fit needs that component to fit the high-energies.