Test differenze fit polyfit vs curve_fit, con e senza errore + minuit al fondo¶

import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as opt

v,xl,dxl = np.loadtxt('planck.dat', unpack=True)

x=3.0/(1.6*v*1000)*1e5


fig, ax = plt.subplots(2,2, figsize=(15,8))

# PAD 0
# Fit
pars0, cov0 = np.polyfit(x, xl, 1, w = None, cov=True)
errors0 = np.sqrt(np.diag(cov0))
# Plot
ax[0][0].set_title('numpy polyfit, no errore')
ax[0][0].plot(x, xl,'o', color='k')
ax[0][0].plot(x, pars0[0]*x + pars0[1] , '--r', label=f'Slope = {pars0[0]:.3f} +- {errors0[0]:.3f}')
ax[0][0].legend(fontsize=14, loc=2)

# Fit
pars1, cov1 = np.polyfit(x, xl, 1, w = 1/dxl, cov='unscaled')
errors1 = np.sqrt(np.diag(cov1))
ax[0][1].set_title('numpy polyfit, pesi = 1/dxl')
ax[0][1].errorbar(x, xl,yerr=dxl, fmt='o', color='k')
ax[0][1].plot(x, pars1[0]*x + pars1[1] , '--r', label=f'Slope = {pars1[0]:.3f} +- {errors1[0]:.3f}')
ax[0][1].legend(fontsize=14, loc=2)

def line(x, m, q) :
    return m*x + q
    
# Fit 
pars2, cov2 = opt.curve_fit(line,x,xl)
errors2 = np.sqrt(np.diag(cov2))
ax[1][0].set_title('scipy.optimize curve_fit, no errori')
ax[1][0].plot(x, xl,'o', color='k')
ax[1][0].plot(x, pars2[0]*x + pars2[1] , '--r', label=f'Slope = {pars2[0]:.3f} +- {errors2[0]:.3f}')
ax[1][0].legend(fontsize=14, loc=2)


# Fit 
pars3, cov3 = opt.curve_fit(line,x,xl, sigma=dxl, absolute_sigma=True)
errors3 = np.sqrt(np.diag(cov3))
ax[1][1].set_title('scipy.optimize curve_fit, errori = dxl')
ax[1][1].errorbar(x, xl,yerr=dxl, fmt='o', color='k')
ax[1][1].plot(x, pars3[0]*x + pars3[1] , '--r', label=f'Slope = {pars3[0]:.3f} +- {errors3[0]:.3f}')
ax[1][1].legend(fontsize=14, loc=2)

plt.show()

Minuit¶

from iminuit import Minuit

def line(x, m, q):
    return m*x + q

def least_squares(m, q):
    return sum((xl - line(x, m, q)) ** 2 / dxl**2)


minuit = Minuit(least_squares, pedantic=False)
minuit.migrad() # finds minimum of least_squares function
minuit.hesse()  # computes errors

for p in minuit.parameters:
    print("{} = {} +/- {}".format(p, minuit.values[p], minuit.errors[p]))

m = 7.022014854577 +/- 1.5914567074169368
q = -1.335400427512674 +/- 9.813945188185663

slope = minuit.values['m']
slope_err = minuit.errors['m']

q = minuit.values['q']
q_err = minuit.errors['q']


plt.title('Minuit')
plt.errorbar(x, xl,yerr=dxl, fmt='o', color='k')
plt.plot(x, slope*x + q , '--r', label=f'Slope = {slope:.3f} +- {slope_err:.3f}')
plt.legend(fontsize=14, loc=2)
plt.show()

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	m	7.0	1.6
1	q	-1.3	9.8