Scipy OptimizeWarning: Covariance of the parameters could not be estimated when trying to fit function to data

I'm trying to plot some data with a non-linear fit using the function:

kB and Bv being a constant while J'' is the independent variable and T is a parameter.

I tried to do this in Python:

#Constants
kb = 1.380649e-23
B = 0.3808

#Fit function
def func(J, T):
    return (B*(2*J+1) / kb*T * np.exp(-B*J*(J+1)/kb*T))
popt, pcov = optimize.curve_fit(func, x, y)
print(popt)

plt.plot(x, y, "o")
plt.plot(x, func(j, popt[0]))

But this results in the warning "OptimizeWarning: Covariance of the parameters could not be estimated warnings.warn('Covariance of the parameters could not be estimated'" and the parameter turns out to be 1.

This is the data:

x = [2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46,48]
y = [0.185,0.375,0.548,0.695,0.849,0.931,0.996,0.992,1.0,0.977,0.927,0.834,0.691,0.68,0.575,0.479,0.421,0.351,0.259,0.208,0.162,0.135,0.093,0.066]

Why?

You have three problems in your model:

• Priority of operations are not respected, as pointed by @jared, you must add parenthesis around the term kb * T;
• Target dataset is normalized to unity where it is its area which is unitary (it is a distribution). therefore such normalization prevents the model to converge. Adding the missing parameter N to the model allows convergence;
• Scale issue and/or units issue, your Boltzmann constant is expressed in J/K but we don't have any clue about the units of B. Recall the term inside the exponential must be dimensionless. Scaling constant to eV/K seems to provide credible temperature.

MCVE

Bellow a complete MCVE fitting your data:

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize, integrate
from sklearn import metrics

#kb = 1.380649e-23    # J/K
kb = 8.617333262e-5   # eV/K
B = 0.3808            # ?

def model(J, T, N):
    return N * (B * (2 * J + 1) / (kb * T) * np.exp(- B * J * (J + 1)/ (kb * T)))

J = np.array([2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46,48])
N = np.array([0.185,0.375,0.548,0.695,0.849,0.931,0.996,0.992,1.0,0.977,0.927,0.834,0.691,0.68,0.575,0.479,0.421,0.351,0.259,0.208,0.162,0.135,0.093,0.066])

popt, pcov = optimize.curve_fit(model, J, N, p0=[1e5, 1])
#(array([2.51661517e+06, 2.75770698e+01]),
# array([[1.17518427e+09, 3.25544771e+03],
#        [3.25544771e+03, 6.04463335e-02]]))

np.sqrt(np.diag(pcov))
# array([3.42809607e+04, 2.45858361e-01])

Nhat = model(J, *popt)
metrics.r2_score(N, Nhat)  # 0.9940231921241076

Jlin = np.linspace(J.min(), J.max(), 200)
Nlin = model(Jlin, *popt)

fig, axe = plt.subplots()
axe.scatter(J, N)
axe.plot(Jlin, Nlin)
axe.set_xlabel("$J''$")
axe.set_ylabel("$N(J'')$")
axe.grid()

Regressed values are about T=2.51e6 ± 0.03e6 and N=2.76e1 ± 0.02e1, both have two significant figures. The fit is reasonable (R^2 = 0.994) and looks like:

If your Boltzmann constant is expressed in J/K, temperature T is about 1.6e+25 which seems a bit too much for such a variable.

Check

We can check the area under the curve is indeed unitary:

Jlin = np.linspace(0, 100, 2000)
Nlin = model(Jlin, *popt)
I = integrate.trapezoid(Nlin / popt[-1], Jlin)  # 0.9999992484190978