Pass ncall and iterate in methods that build profile likelihoods
#1008
Comments
|
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 |
|
It is a good idea, thanks. Will do. |
|
I can write up a PR, I don't know how much customizability you'd want though. I think |
|
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. |
|
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 |
|
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 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 |
|
I got 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 With the hack, I'm just getting the hack warning but then the minima duing the scan are valid. |
|
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. |
|
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 |
|
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. |
|
When I try your example (thanks), it fails with an error in the step |
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. |
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. |
|
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. |
|
Yup, I'm on arm64 macos. |
|
Anyway, I get the warnings, so I can try to debug it, thanks for now. |
|
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. |
|
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. |
|
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. |
|
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. |
#1010 there you go |
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. |
|
The C++ code generally does not care a lot about internal consistency. |
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. |
|
You got it, that is the issue. |
|
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 |
|
Thank you Hans for the PR and the tutorial!!. Should be helpful to people struggling with such fits. |
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. |
Hi @HDembinski,
Would you be open to allow
migradoptional arguments as arguments to method that return profile likelihoods likemnprofile. 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.
The text was updated successfully, but these errors were encountered: