[Python] iminuit not converging results

I am using Iminuit to find the sensitivities of parameters. I have signal and background which depends to energy. Before launching IMINUIT, we have to define the energy range and the number binning. We have noticed that more we increase the number of bin (~50000), more upper limit decrease..It seems it not converge. How it is possible? depends on iminuit? how I understand which is the exact value of number binning?

Table is a table with three columns: mass, energy, flux..for each mass there are ~800values of energies

class LLF:

def __init__(self, model, Energy, non, noff, alpha, Tobs):
self.model = model
self.Energy = Energy
self.non = np.asanyarray(non, dtype=np.float64)
self.noff = np.asanyarray(noff, dtype=np.float64)
self.alpha = alpha
self.Tobs = Tobs

def __call__(self, *par):

stat = likelihood(self.Energy, self.non, self.noff, self.alpha, self.Tobs, *par)

return 2 * stat

class BetterLLF(LLF):

def __init__(self, model, Energy, non, noff, alpha, Tobs):
super().__init__(model, Energy, non, noff, alpha, Tobs)
self.func_code = make_func_code(describe(model)[5:])

def likelihood(Energy, non, noff, alpha, Tobs, norm, tilt, log_g, ma):
g = np.power(10,log_g)

nbgd = expectedBkg(Energy, Tobs) * norm * np.power(Energy,tilt)
ngam = expectedSignal(Energy, Tobs, g)
nonpred = ngam + alpha * nbgd

term1 = nonpred + nbgd
with np.errstate(divide="ignore", invalid="ignore"):
term2_ = -non * np.log(nonpred)
condition = non == 0
term2 = np.where(condition, 0, term2_)

with np.errstate(divide="ignore", invalid="ignore"):
term3_ = - noff * np.log(nbgd)
condition = noff == 0
term3 = np.where(condition, 0, term3_)

value = np.sum(term1) + np.sum(term2) + np.sum(term3)

return value


NRUNS=1

RATIO = 2.71

Emin=200

Emax=10000

Ea = np.linspace(Emin, Emax, 1000) ##HERE THE NUMBER OF BINS

alpha = 1
Tobs = 50000 #observation time

# Fitting parameter:

g = 1e-26

log_g = np.log10(g)

columns = ['model', 'g_95','g_test', 'ma']

table = pd.DataFrame(columns=columns)

np.random.seed(int(time.time()))

g_95 = []
g_test = []
model = []
lgma_array = []

for lgma in np.unique(Table['lgma']):

mask = (Table['lgma']==lgma)
Spectrum = interp1d(Table[mask]['lgE'],Table[mask]['Flux'], kind='previous', fill_value='extrapolate')

ma = 10**lgma
g_95_array = []
for run in range(NRUNS):

noff = expectedBkg(Ea, Tobs)
n_model = expectedSignal(Ea, Tobs, g)

non = n_model + noff

llf = BetterLLF(likelihood, Ea, non, noff, alpha, Tobs)
m_test = Minuit(llf, norm=1, tilt=0, log_g=log_g, ma=ma)
m_best = Minuit(llf, norm=1, tilt=0, log_g=log_g, ma=ma)

m_test.errordef=Minuit.LEAST_SQUARES
#m_test.limits['tilt'] = (-1,1)
#m_test.limits['norm'] = (0.5, 1.5)

m_test.fixed["ma"] = True
m_test.fixed["log_g"] = True
R_test = m_test.migrad()


m_best.errordef=Minuit.LEAST_SQUARES
#m_best.limits['tilt'] = (-1,1)
#m_best.limits['norm'] = (0.5, 1.5)

m_best.fixed["ma"] = True
m_best.fixed["log_g"] = True
R1_best = m_best.migrad()

m_best.fixed["log_g"] = False
#m_best.limits['log_g'] = (-28,-21)
R_best = m_best.migrad()

print(R_test.valid )
print(R_best.valid )

if R_best.valid == False:
continue


print(10**m_test.values['log_g'], 10**m_best.values['log_g'])

scale_min = R_best.values["log_g"]
scale_max = scale_min + 0.4

scan_par, scan_value, _ = R_best.mnprofile(vname="log_g", size=200,
bound=(scale_min,scale_max))
scan_value = scan_value-R_test.fmin.fval



sigma = 0.6
max_iter = 100
counter = 0
while scan_value[-1]-RATIO <=0 and counter<max_iter:
scale_max = scale_min + sigma
scan_par, scan_value, _ = R_best.mnprofile(vname="log_g", size=200,
bound=(scale_min,scale_max))
scan_value = scan_value-R_test.fmin.fval
sigma = sigma + 0.2
counter = counter + 1

xvals = scan_par
yvals = scan_value - RATIO

scale_found = brentq(
interp1d(xvals, yvals, kind="quadratic"),
scale_min, scale_max,
maxiter=100,
rtol=1e-5,
)


print(run, counter, ma, np.power(10,scale_found))
g_95_array.append(np.power(10,scale_found))


g_95.append(np.quantile(g_95_array, [0.05,0.32,0.5,0.68,0.95])[2])
g_test.append(m_test.values['log_g'])
model.append(0)
lgma_array.append(lgma)

table = pd.DataFrame({'model':model, 'g_95':g_95, 'g_test':g_test, 'lgma':lgma_array})

Continue reading...