I’m trying to fit modeling data as a curve using Scipy’s optimize curve_fit to some scattered data, but I got a negative value of (bX) parameter that makes no sense. So I’m not sure what is wrong with the codes.

# Define linear quadratic curve
def linq(x, aparam, bparam):
    return np.exp(-(aparam*x+ bparam*x**2))
#new model
NX=46
def abX(x, aXparam, bXparam):
    return np.exp(-(aXparam + bXparam)*x)*(1 + bXparam*x/NX)**NX
# define result array
avalues =[]
davalues =[]
bvalues=[]
dbvalues =[]
Chi2=[]

# Model parameters
Dmax = 10.
a1=0.740; b1=.1100
# ensemble size
Nensemble = 10

for iexp  in range(1, Nensemble+1):
# Create the artificial dataset
    nobs = int(Dmax + .5)
    x = np.arange(.5,nobs)
    N = linq(x,a1,b1)
    Nfluct = 0.3*N*np.random.normal(size=nobs)
    N = N + Nfluct
    sig = np.zeros(nobs) + 0.3*N
# Fit the curve
    start = (0,1.)
    popt, pcov = curve_fit(abX,x,N,sigma = sig,p0 = start,absolute_sigma=True)
  
# add result to result vectors
    avalues=avalues+[popt[0]]
    bvalues=bvalues+[popt[1]]
    davalues=davalues+[np.sqrt(pcov[0,0])]
    dbvalues=dbvalues+[np.sqrt(pcov[1,1])]

print("mean aX value =",np.mean(avalues),"+/-",np.std(avalues)/np.sqrt(Nensemble-1),", gen =",a1) 
print("mean bX value =",np.mean(bvalues),"+/-",np.std(bvalues)/np.sqrt(Nensemble-1),", gen =",b1)

mean aX value = 0.9829722097543518 +/- 0.019840535134131847 , gen = 0.74
mean bX value = -2.3048945875855527 +/- 0.024140253716692855 , gen = 0.11