CODEWITHSUNDEEP
pythoncurve-fittingscipy-optimizeleast-squaresnon-linear-regression
yusra zabarmawi
yusra zabarmawi

Reputation: 1

Why do I get negative values for parameters in a curve fitting?

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

Upvotes: 0

Views: 104

Answers (1)

jlandercy
jlandercy

Reputation: 11042

TL; DR

There are multiple minima for the second model, changing the initial parameter for (0., 5.) allow the algorithm to reach the second minima where b>0.

Or as suggested by reinderien add bounds to filter out the first minima.

Checking model compliance

First lets check your second model can fit the first one:

a1 = 0.740
b1 = 0.1100

def model1(x, a, b):
    return np.exp(-(a * x + b * x **2))

def model2(x, a, b, N=46):
    return np.exp(- (a + b) * x) * (1 + b * x / N) ** N

xlin = np.linspace(0.5, 9.5, 200)
ylin = model1(xlin, a1, b1)

popt, pcov = optimize.curve_fit(model2, xlin, ylin, p0=[0., 1.])
# (array([ 0.7524258 , -2.93622762]),
#  array([[2.51911174e-07, 1.96469005e-06],
#         [1.96469005e-06, 1.92412517e-05]]))

We can confirm that we have correct fitness between both model.

enter image description here

Checking optimality

Now lets check this is the optimal set of parameters:

def loss_factory(model, x, y):
    @np.vectorize
    def wrapper(a, b):
        return 0.5 * np.sum(np.power(y - model(x, a, b), 2))
    return wrapper

alin = np.linspace(0, 2, 300)
blin = np.linspace(-5, 5, 300)
A, B = np.meshgrid(alin, blin)

loss = loss_factory(model2, xlin, ylin)

SSE = loss(A, B)

enter image description here

We can see there are two minima in the explored parameters region. Where the lower one (negative b) is found because of the initial parameters.

To reach the second minima, just change initial parameters:

popt, pcov = optimize.curve_fit(model2, xlin, ylin, p0=[0., 5.])
# (array([0.72716386, 3.44478853]),
#  array([[ 2.99618275e-07, -2.70844576e-06],
#         [-2.70844576e-06,  2.85220271e-05]]))

Or add bounds:

popt, pcov = optimize.curve_fit(model2, xlin, ylin, p0=[0., 1.], bounds=(0., np.inf))
# (array([0.72716384, 3.44478875]),
#  array([[ 2.99617959e-07, -2.70844633e-06],
#         [-2.70844633e-06,  2.85220715e-05]]))

The very nature of error landscape SSE is due to the dataset (here strictly following the model2) and the existence of multiple minima is most probably an effect of the even power (N=46).

Upvotes: 1

Related Questions