pythonnumpyscipycurve-fittinghyperbolic

Fitting Hyperbolic Cosine curve in Python

Now I want to fit in one bump of hyperbolic cosine curve into the following X and Y data:

xData = np.array([1.7, 8.8, 15, 25, 35, 45, 54.8, 60, 64.7, 70])
yData = np.array([30, 20, 13.2, 6.2, 3.9, 5.2, 10, 14.8, 20, 27.5])

Here's what I have done so far but I am not getting the expected result and I have no idea what I am doing wrong:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import scipy.interpolate as inp

xData = np.array([1.7, 8.8, 15, 25, 35, 45, 54.8, 60, 64.7, 70])
yData = np.array([30, 20, 13.2, 6.2, 3.9, 5.2, 10, 14.8, 20, 27.5])

def model_hcosine(x, a, b, c):
    return a * np.cosh(x/b) + c

poptcosh, pcovcosh = curve_fit(model_hcosine, xData, yData, p0=[min(yData), max(xData), max(yData)])

aapopt, bbopt, cccopt = poptcosh
xCoshModel = np.linspace(min(xData), max(xData), 100)
yCoshModel = model_hcosine(xCoshModel, aapopt, bbopt, cccopt)

plt.scatter(xData, yData)
plt.plot(xCoshModel, yCoshModel, 'b-')

plt.show()
Solution

  • @WarrenWeckesser is correct, you need to account for the translation within the cosh function. You can add an additional parameter d to the model, and give it an initial condition of 0 in the optimizer. Then you unpack the optimal coefficients and plug them into the model before plotting. I got the following

    import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import scipy.interpolate as inp

xData = np.array([1.7, 8.8, 15, 25, 35, 45, 54.8, 60, 64.7, 70])
yData = np.array([30, 20, 13.2, 6.2, 3.9, 5.2, 10, 14.8, 20, 27.5])

def model_hcosine(x, a, b, c, d):
    return a * np.cosh((x-d)/b) + c

poptcosh, pcovcosh = curve_fit(model_hcosine, xData, yData, p0=[min(yData), max(xData), max(yData), 0])

aapopt, bbopt, cccopt, ddopt = poptcosh
xCoshModel = np.linspace(min(xData), max(xData), 100)
yCoshModel = model_hcosine(xCoshModel, aapopt, bbopt, cccopt, ddopt)

plt.scatter(xData, yData)
plt.plot(xCoshModel, yCoshModel, 'b-')

plt.show()

    result