Search code examples
pythonscipy-optimize

Bad curve_fit result

I'm trying to fit data with the next custom function:

But when plotting iC and iD with the regression-adjusted constants, the results I get are as follows.

enter image description here

Where RS, CDL, B and E are the fitting constants from the regression.

My data is taken from a CSV with two columns, where X represents time and Y represents iT.

Here are some of the dataframe data:

               X         Y
903   1147.400171  0.004707
904   1148.400371  0.000036
905   1149.400371  0.110646
906   1150.400571  0.000008
907   1151.400571  0.034051
...           ...       ...
1199  1443.400764 -0.001796
1200  1444.400764 -0.001924
1201  1445.400763 -0.001779
1202  1446.400763 -0.001818
1203  1447.399963 -0.001818

My code to call the function looks as follows:

def model_f(x, E, Rs, CDL, B):
    d = (math.pow(E/Rs, -(x/Rs*CDL))) + (B/math.sqrt(x)) + E
    return d

model_f = np.vectorize(model_f)

def get_difusive(B, t):
        return B/math.sqrt(t)

def get_capacitive(t, E,Rs,CDL):
        return math.pow((E/Rs), -(t/Rs*CDL))

                
x_data = df['X'].values
y_data = df['Y'].values

popt, pcov = curve_fit(f=model_f,xdata= x_data, ydata= y_data, p0=(1,1,1,1), bounds=(1,3))
E, Rs, CDL, B = popt

y = np.empty(len(x_data))
for i in range(len(x_data)):
    y[i] = get_difusive(B, x_data[i])
    
plt.plot(x_data, y, 'r--', label='Difusiva')

y = np.empty(len(x_data))
for i in range(len(x_data)):
    y[i] = get_capacitive(x_data[i], E, Rs, CDL)
    
plt.plot(x_data, y, 'b--', label='Capacitiva')

    
plt.plot(x_data,y_data, label='Original data')
plt.legend()
plt.show()

It is clear that the model does not fit the data correctly. I have tried to change the settings that curve_fit gives me and I have also tried doing a non-linear regression with other models like lmfit, but it hasn't worked.

Solution

  • It's not a perfect fit, but it's a fit. Your call to pow() was incorrect, and you need bounds. You also need to filter away your bad data points, and subtract the initial time offset.

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


def diffusive(t: np.ndarray, t0: float, E: float, Rs: float, Cdl: float, B: float) -> np.ndarray:
    return B/np.sqrt(t - t0)  # + E


def capacitive(t: np.ndarray, t0: float, E: float, Rs: float, Cdl: float, B: float) -> np.ndarray:
    return E*(1 + np.exp(-(t - t0)/Rs/Cdl)/Rs)


def decay_model(t: np.ndarray, t0: float, E: float, Rs: float, Cdl: float, B: float) -> np.ndarray:
    i = (capacitive(t, t0, E, Rs, Cdl, B)
        + diffusive(t, t0, E, Rs, Cdl, B))
    return i


df = pd.read_csv('raw_data.csv', names=['time', 'current'], skiprows=1)
df.time -= df.time.min()
df = df[
    (df.time > 25) |
    (df.current > 2.5e-3)
].loc[1:, :]
popt, *_ = curve_fit(
    f=decay_model, xdata=df.time, ydata=df.current,
    p0=(0.1,  0.1, 1e-2, 1e2, 0),
    bounds=(
        ( 0, 1e-6, 1e-3,   1,   0),
        (10,    1,    1, 1e3,   1),
    ),
)
print(popt)

fig, ax = plt.subplots()
plt.plot(df.time, df.current, label='Original data')
ax.semilogx(df.time, diffusive(df.time, *popt), label='Diffusive fit')
ax.semilogx(df.time, capacitive(df.time, *popt), label='Capacitive fit')
ax.semilogx(df.time, decay_model(df.time, *popt), label='Total fit')

plt.legend()
plt.show()

    [2.00019972e+00 1.39365575e-03 4.26528691e-02 1.20252864e+02
 3.67268296e-05]

    fit