Anyway to fit this model to these data points?

Question

I have been tryin to fit these data points to the exponential model y=ae^(px+qx^2). The reason I need this model fit is due to the paper "Temperature dependence of the functional response" by Göran Englund 2011 where they fit this model on data points as described in one of the figures. I need to fit this with some data I found doing literature review for an ODE model I am working on in my dissertation. However no matter what I do, I cannot fit it using the curvefit() function, even when my guess p0 is very close. Below is my code and a figure of my fit when I just guessed at the fitted parameters a,p and q until I got a super close fit.

k = 0.00008617
def NormalizeData(data):
    return (data - np.min(data)) / (np.max(data) - np.min(data))

def C_to_A(C):
    return -1/(k*(C+273))

def model(x, a, p,q):
    return a*np.exp(p*x+q*x**2)

# Input data
x = C_to_A(np.array([11, 11, 11,15,15,15,20,20,20,25,25,25,29,29,29]))
y = NormalizeData(np.array([2.048, 1.56, 1.18, 2.6116,2.35,2.1036,2.97, 2.97, 2.97, 2.463,2.05,1.6679,1.825,1.0939,0.534]))
# do the fit with some initial values
popt, pcov = curve_fit(model, x, y, p0=(9.7*10**-139.95, -16, -0.2))

# prepare some data for a plot
xx = np.linspace(-60, 0)
yy = model(xx, *popt)

plt.plot(x, y, 'o', xx, yy)
plt.title('Exponential Fit')
plt.ylim(0,1.1)
plt.xlim(-50,32)
plt.grid()
plt.show()
print(popt)

If I just plot the model with the data points and my guessed coefficients a=9.7*10**-139.95, p=-16, q=-0.2 I get the picture below

Which is super close to what I want. Any suggestions on how I may approach this fit using the model I mentioned? Anything will help!

Answer 1

y=ae^(px+qx^2) implies log(y) = c + px + qx^2, where c=log(a),p,q are parameters. You can estimate the parameters using linear regression of y on x and x^2.

import numpy as np
import pandas as pd
from statsmodels.api import OLS
import matplotlib.pyplot as plt


def C_to_A(C, k=0.00008617):
    return -1/(k*(C+273))

# choice of k is arbitrary, change as you wish
x = C_to_A(np.array([11, 11, 11,15,15,15,20,20,20,25,25,25,29,29,29]), k=1)
y = np.array([2.048, 1.56, 1.18, 2.6116,2.35,2.1036,2.97, 2.97, 2.97, 2.463,2.05,1.6679,1.825,1.0939,0.534])


# want to regress log(y) on x, x^2, constant
# create data for this
data = pd.DataFrame({"y":y, "x":x})
data["log_y"] = np.log(data.y)
data["x2"] = data.x**2
data["const"] = 1
data

Estimate the parameters using OLS from statsmodels:

Y = data.log_y
X = data[["x", "x2", "const"]]

ols = OLS(Y, X, hasconst=True).fit()
ols.summary()

For plotting, ols's fitted values property is handy. We convert log(y) to y using exp:

plt.scatter(x, y)
plt.plot(x, np.exp(ols.fittedvalues))