scipy.optimize.curve_fit giving same values as that of initial values

Question

I am trying to fit a curve to the data that I have. Seems like an easy thing to do but whenever I fit the curve, the optimization doesn't seem to work. The curve_fit optimization just spits out whatever initial values I give. If I don't give any initial values for the parameters, it just gives [1, 1, 1] which is the default initial guess used by curve_fit. Can anyone point out the mistake that I am making?

Thanks,

Following is the code and data that I am using:

import pandas as pd
import numpy as np
from scipy.optimize import curve_fit

x = [0.3, 0.6, 0.8, 1, 1.1, 1.5, 1.9, 2.4, 3.2, 12.6]
y = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100]

df = pd.DataFrame({"X": x, "Y": y})

# Function
def sb_model(x, xmax, xc, b_val):
    return (1/(1+((np.log(xmax/x)/np.log(xmax/xc))**b_val)))*100

popt, pcov = curve_fit(f=sb_model, xdata=df["X"], ydata=df["Y"], p0 = [10.0, 0.9, 3.2])
print(popt)

Output:

[10.   0.9  3.2]

Answer 1

I did manage to get your model to fit your data. I believe you were trying to fit it without looking at the resulting graphs. That's the first thing to do when doing fits. Often in non-linear fitting, you have to provide a set of initial points as close to the result as possible, so the algorithm has an easier time reaching it.

Doing just that, I noticed your model was much too low compared to the data and I couldn't get it to "go up" by modifying any parameter. I ended up adding a "scale" parameter. Next, I noticed you were getting errors because there were negative values in your logs. I fixed that by setting the minimum value of your parameters to zero. I ended up doing this first in lmfit, then tested with curve_fit, and both worked. Here's the code and the graphs.

import pandas as pd
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

x = [0.3, 0.6, 0.8, 1, 1.1, 1.5, 1.9, 2.4, 3.2, 12.6]
y = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100]

df = pd.DataFrame({"X": x, "Y": y})

# Function
def sb_model(x, xmax, xc, b_val, scale):
    return scale/(1+((np.log(xmax/x)/np.log(xmax/xc))**b_val))

plt.plot(x, y)
p0 = np.array([100, 0.1, 1, 100])
model = sb_model(x, *p0)
plt.plot(x, model)

Gave me this model and data, which are relatively close:

I fitted it with

popt, pcov = curve_fit(f=sb_model, xdata=df["X"], ydata=df["Y"], p0 = p0, maxfev=800, bounds=(0, 1E5))

And got this:

popt was array([271.35537164, 1.21550569, 10.44652752, 100.36436418])

If you're interested in the lmfit code, here's the code. Its output was essentially the same, but it also provides a standard error (which, BTW, is 500% for xmax and 90% for b_val):

from lmfit import Parameters, minimize
p = Parameters()
p.add('xmax', value=100, min=0)
p.add('xc', value=1, min=0)
p.add('b_val', value=1)
p.add('scale', value=100, vary=True)

def residual(pars, x, data):
    xmax = pars['xmax']
    xc = pars['xc']
    b_val = pars['b_val']
    scale = pars['scale']
    model = sb_model(x, xmax, xc, b_val, scale)
    return model - data
  
out = minimize(residual, p, args=(x, y))

Adding the constraints requested:

from lmfit import Parameters, minimize

# xmax > xc
# xmax - xc > 0
# xmax - xc = delta; delta > 0
# xc = xmax - delta

p = Parameters()
acceptable_min = 1e-10 # Just so the values don't reach 0
p.add("xmax", value=100, min=acceptable_min, max=max(x))
p.add("delta", value=110, min=acceptable_min)
p.add("xc", min=acceptable_min, value=10, expr="delta - xmax")
p.add("b_val", value=1, min=acceptable_min)

def sb_model(x, xmax, xc, b_val):
    first_log = np.log(xmax / x)
    second_log = np.log(xmax / xc)              
    division = (first_log / second_log)
    # The suggestion by @Reinderen fixed the problem I was having here
    power =  np.sign(division) * np.abs(division) ** b_val
    divisor = 1 + power
    val = 100 / divisor # Added the '100' factor here directly
    
    # I was having some problems with nans so I added this to monitor them.
    #print(f'{xmax.value=}, {xc.value=}, {b_val.value=}')
    #print(f'{first_log=}')
    #print(f'{second_log=}')
    #print(f'{division=}')
    #print(f'{power=}')
    #print(f'{divisor=}')
    #print(f'{val=}')
    #print()
    return val


def residual(pars, x, data):
    xmax = pars["xmax"]
    xc = pars["xc"]
    b_val = pars["b_val"]
    model = sb_model(x, xmax, xc, b_val)
    return model - data

out = minimize(residual, p, args=(x, y))

plt.plot(x, y)
plt.plot(x, y + out.residual)

If I remove the maximum value of xmax, I get: