我正在尝试使用 minuit 拟合曲线，但似乎数据未正确拟合曲线。我怀疑它与函数中的日志有关，因为一旦我从 chisquare 函数中删除 chiaf，值就会正确拟合。问题还在于值 am18 为负值，它对应于一个质量，并且由于所有其他值都是正值，所以这个值不是正值真的很奇怪。
import numpy as np
from iminuit import Minuit
#some constants
hc = 197.3269788
mssph = 685.78
#some data and errors to fit
muds = np.array([
81.96100485, 64.96427609, 38.15520137, 45.75992993,
27.38344029, 23.41742996, 18.73586921, 18.07486749,
9.20589292, 4.83878931, 72.17070899, 71.08083681,
39.57647386, 31.63373626, 30.65363538, 26.50584842,
19.33150992, 13.45482499, 3.86943831, 68.76598171,
66.7933219 , 38.86139831, 38.61246743, 38.6500419 ,
16.62670414, 16.24156579, 9.16347451, 52.48804692,
51.95712296, 23.86333781, 23.81815279, 12.44282442
]) #xdata
Fr = np.array([
127.31527434, 122.72790265, 110.26449558, 112.75717699,
104.81830088, 104.35746903, 101.32016814, 100.54513274,
96.87942478, 92.98330088, 124.9736053 , 122.52414305,
112.47114172, 108.74591788, 107.34258013, 108.00597616,
102.18850331, 100.04522384, 91.47210596, 128.18641113,
122.15516847, 108.23229985, 109.85263369, 107.69218856,
99.14042658, 98.0902102 , 99.1104204 , 112.47678261,
110.39126087, 98.373 , 98.97391304, 95.01495652
])#ydata
Zsrgi = np.array([1.47, 1.5, 1.54, 1.58])
dzs = np.array([ 3.60555128, 3.60555128, 4.24264069, 1.41421356])
Za = np.array([ 0.9468, 0.9632, 0.9707, 0.9756])
Zaey = np.array([56.0, 53.0, 35.0, 15.0])
dmfel = np.array([
2743.82452989, 1388.6310737 , 689.77566743, 652.7213305 ,
348.31861265, 309.31427565, 134.12993165, 157.95193019,
60.83965035, 19.80905202, 3482.5565895 , 3378.1685061 ,
949.07161448, 570.68604168, 492.42222648, 801.32561084,
362.30696375, 152.61749342, 26.39235787, 2733.32712033,
2329.20044751, 788.45797536, 754.06472854, 1098.9552043 ,
279.63806125, 213.46692781, 169.87677303, 2181.08350448,
2018.45039786, 696.73902333, 763.51325268, 416.74304481
])
af = np.array([
0.0575465 , 0.05547301, 0.04983955, 0.05096624, 0.04737787,
0.04716958, 0.04579672, 0.0454464 , 0.0437895 , 0.04202845,
0.04717754, 0.04625286, 0.04245786, 0.04105158, 0.04052182,
0.04077226, 0.03857616, 0.03776707, 0.03453072, 0.0414683 ,
0.0395172 , 0.03501315, 0.03553733, 0.03483842, 0.03207193,
0.03173218, 0.03206222, 0.03104359, 0.03046799, 0.02715095,
0.0273168 , 0.02622413
])
afe = np.array([
0.00039571, 0.00046823, 0.00034898, 0.00071523, 0.00052648,
0.00051809, 0.00072373, 0.00059257, 0.00054912, 0.00065308,
0.000283 , 0.00028516, 0.00031909, 0.0003537 , 0.00035294,
0.0002573 , 0.00034244, 0.00034139, 0.00055726, 0.00038002,
0.00044435, 0.00050214, 0.00052077, 0.00037143, 0.00036278,
0.00045654, 0.00039068, 0.00025011, 0.00024068, 0.00031897,
0.00037707, 0.0003581
])
amssr = np.array([
0.35006, 0.34592, 0.33967, 0.3175 , 0.31341, 0.34965,
0.33393, 0.31033, 0.30807, 0.32813, 0.29895, 0.26546,
0.29559, 0.29298, 0.26026, 0.29264, 0.29094, 0.29074,
0.25928, 0.3888 , 0.23213, 0.23284, 0.22768, 0.20825,
0.22591, 0.20437, 0.22353, 0.2036 , 0.18987, 0.20088,
0.18743, 0.18801
])
amudr = np.array([
0.01737619, 0.01377279, 0.00808912, 0.00970136, 0.00580544,
0.00496463, 0.00397211, 0.00383197, 0.0019517 , 0.00102585,
0.01227267, 0.01208733, 0.00673 , 0.00537933, 0.00521267,
0.00450733, 0.00328733, 0.002288 , 0.000658 , 0.00950714,
0.00923442, 0.00537273, 0.00533831, 0.00534351, 0.0022987 ,
0.00224545, 0.00126688, 0.00588165, 0.00582215, 0.00267405,
0.00266899, 0.0013943
])
amudre = np.array([
8.75479519e-04, 7.34296720e-04, 4.51927448e-04,
5.98343236e-04, 3.81187530e-04, 3.18993055e-04,
4.05218892e-04, 3.28835732e-04, 2.14302222e-04,
1.80198850e-04, 6.27776371e-04, 6.18193827e-04,
3.56742512e-04, 3.04857443e-04, 3.08046099e-04,
2.32852126e-04, 1.84120801e-04, 1.62566426e-04,
7.02801597e-05, 5.87470469e-04, 5.70265978e-04,
3.42904960e-04, 3.42757519e-04, 3.31231953e-04,
1.52175013e-04, 1.61067797e-04, 8.21035930e-05,
1.40038557e-04, 1.42199518e-04, 7.54950814e-05,
7.13899224e-05, 3.64415759e-05
])
mse = np.array([
0.00049366, 0.00064125, 0.00045706, 0.00094021, 0.0009005 ,
0.00074027, 0.00093984, 0.00068264, 0.00086371, 0.00083934,
0.0004669 , 0.00043174, 0.00057567, 0.00059908, 0.00064938,
0.0005022 , 0.00082 , 0.00090521, 0.00106902, 0.00073246,
0.00075007, 0.00098326, 0.0011934 , 0.0006229 , 0.00081216,
0.00058873, 0.00052038, 0.00055326, 0.0005316 , 0.00088051,
0.00105233, 0.00053235
])
a = np.array([ 0.0904, 0.0755, 0.0647, 0.0552])
B = 2.61
Fc = 92.8
def Ffitr(X, s, pp, k, B=2.61, Fc=92.8, mu=770, Za=0.9468): #fit function
logBXmu = np.log( (2 * B * X) / mu**2)
temp1 = (2 * B * X) / (4 * pi * Fc)**2
temp2 = temp1 * (afij[0] +
afij[1] * logBXmu)
temp3 = temp1**2*( afij[2] +
afij[3] * logBXmu +
afij[4] * logBXmu**2)
temp4 = temp1**3 * (afij[5] +
afij[6] * logBXmu +
afij[7] * logBXmu**2 +
afij[8] * logBXmu**3)
return Fc * pp * (1 + k * s) * (1 + temp2 + temp3 + temp4)
#chisquare
def chiwork( a0, a1, a2, a3,
b0, b1, b2, b3, b4, b5, b6, b7, b8, b9,
b10, b11, b12, b13, b14, b15, b16, b17, b18, b19,
b20, b21, b22, b23, b24, b25, b26, b27, b28, b29,
b30, b31,
z0, z1, z2, z3, za0, za1, za2, za3,
am0, am1, am2, am3, am4, am5, am6, am7, am8, am9,
am10,am11, am12, am13, am14, am15, am16, am17, am18, am19,
am20, am21, am22, am23, am24, am25, am26, am27, am28, am29,
am30, am31,
k, y1):
B = 2.61
Fc = 92.8
hc = 197.3269788
mu = 770
ana = np.array([a0, a1, a2, a3])
zs = np.array([z0, z1, z2, z3])
za = np.array([za0, za1, za2, za3])
am = np.array([am0, am1, am2, am3, am4, am5, am6, am7, am8, am9,
am10, am11, am12, am13, am14, am15, am16, am17, am18, am19,
am20, am21, am22, am23, am24, am25, am26, am27, am28, am29,
am30, am31])
#am[numpy.isnan(am)] = 10**(-50)
am[am <= 0] = 10**(-50)
bss = np.array([b0, b1, b2, b3, b4, b5, b6, b7, b8, b9,
b10, b11, b12, b13, b14, b15, b16, b17,b18,b19,
b20, b21, b22, b23, b24, b25, b26, b27, b28, b29,
b30, b31])
betalat = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1, 1, 1, 1, 1, 1, 1, 1, 1, 2,
2, 2, 2, 2, 2, 2, 2, 3, 3, 3,
3, 3]
chia = sum( ( (ana - a) / (ast * 10**(-4) ) )**2) +\
sum( ( (bss * a[betalat] / hc - amssr) / mse )**2)
chiz = sum( ( (zs - 1 / Zsrgi) / (dzs * 10**(-2) ) )**2) +\
sum( ( (za - Za) / (Zaey * 10**(-4) ) )**2)
chiamd = sum( ( (amudr - (1 + y1 * ana[betalat]**2) *
ana[betalat] * am * zs[betalat] / hc) / amudre)**2)
chiaf = sum( ( ( af - Ffitr(am, bss**2 - mssph**2, ana[betalat] / 0.95, k) / hc) / afe)**2)
return chia + chiz + chiamd + chiaf
#minuit part
mc=minuit.Minuit(
chiwork, a0=0.09, z0=0.68, za0=0.95, am0=62,
b0=755, k=5.124105506705958e-07, y1=1,
error_a0=0.02, error_z0=0.02,
error_za0=0.02, error_am0=13, error_b0=4,
error_k=10**(-8), error_y1=0.1,errordef=1
)
fminc, paramc = mc.migrad()
print(mc.values)
#gtting values to plot
asp = []
bs = []
Zs = []
zas = []
mdus = []
for i in range(4):
asp.append(mc.values[i])
for i in range(4, 36):
bs.append(mc.values[i])
for i in range(36, 40):
Zs.append(mc.values[i])
for i in range(40, 44):
zas.append(mc.values[i])
for i in range(44,76):
mdus.append(mc.values[i])
anna = np.array(asp)
bna = np.array(bs)
Zsn = np.array(Zs)
zasn = np.array(zas)
mduss = np.array(mdus)
mnew = np.linspace(min(mduss), max(mduss), 100)
bnew = np.linspace(min(bna), max(bna), 100)
annas = np.linspace(min(anna), max(anna), 100)
zasns = np.linspace(min(zasn), max(zasn), 100)
zsn = np.linspace(min(Zsn), max(Zsn), 100)
plt.figure(figsize=( 14, 8.5))
plt.ylim(90, 200)
plt.plot(muds, Fr, "b^")
plt.plot()
plt.plot(mnew, Ffitr(mnew, bnew**2 - mssph**2, annas / zasns, mc.values["k"]) / annas)
我得到什么：
