Use numpy.curve_fit to on a function with summation terms

Viewed 136

I have a set of x and y data and I am trying to fit them to an equation of the form:

enter image description here

where a and b are the parameters I am trying to solve for.

I defined my fitting function as follow:

def model(x,a,b): x = x[...,np.newaxis] f = np.sum(a*(x**2*b**2)/(1+x**2*b**2), axis = -1) return f

but since a and b are arrays of values I am not sure how to use numpy.curve_fit

This is what I wrote so far:

popt,pcov = curve_fit(lambda x,*params: model(x,a0,b0),
                  x_data,y_data)

Where a0 and b0 are guess 1-D arrays containing 19 elements (same dimensions of my data point set).

The code is only outputting a single 1-D array of 19 elements where I need 2 1-D arrays, one for a and b.

Any suggestions?

x_data = np.array([2.46e-3,4.59e-3,7.46e-3,
              1.23e-2,2.20e-2,3.38e-2,7.76e-2,
              1.33e-1,2.78e-1,6.74e-1,
              1.44e0,3.40e0,8.14e0,
              1.72e1,3.94e1,8.68e1,
              2.55e2,7.62e2,2.03e3])

y_data = np.array([1.18e1,3.70e1,7.13e1,
            1.30e2,2.61e2,4.19e2,9.14e2,
            1.55e3,2.91e3,5.36e3,8.60e3,
            1.40e4,2.28e4,3.32e4,4.69e4,6.46e4,9.52e4,
            1.35e5,1.73e5])

def model(x,a,b):
  x = x[...,np.newaxis]
  f = np.sum(a*(x**2*b**2)/(1+x**2*b**2), axis = -1)
  return f



a0 = np.ones(19)
b0 = np.ones(19)



popt,pcov = curve_fit(lambda x,*params: fit(x,a0,b0),x_data,y_data)

print(popt)
1 Answers

There's probably room for improvement in this solution but you might want to consider using scipy.optimize.minimize for your solution:

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

x_data = np.array([2.46e-3,4.59e-3,7.46e-3,
          1.23e-2,2.20e-2,3.38e-2,7.76e-2,
          1.33e-1,2.78e-1,6.74e-1,
          1.44e0,3.40e0,8.14e0,
          1.72e1,3.94e1,8.68e1,
          2.55e2,7.62e2,2.03e3])

y_data = np.array([1.18e1,3.70e1,7.13e1,
        1.30e2,2.61e2,4.19e2,9.14e2,
        1.55e3,2.91e3,5.36e3,8.60e3,
        1.40e4,2.28e4,3.32e4,4.69e4,6.46e4,9.52e4,
        1.35e5,1.73e5])

def model(x,a,b):
    x = x[...,np.newaxis]
    f = np.sum(a*np.exp(-x/b),axis = -1)
    return f

a0 = np.ones(19)
b0 = np.ones(19)
coefs = np.vstack([a0,b0])

def objectiveFunction(coefs,x,y):
    a = coefs[:19]
    b = coefs[19:]
    modeled = model(x_data,a,b)
    sse = np.sum((y_data-modeled)**2)
    return sse

solution = minimize(objectiveFunction,coefs,args=(x_data,y_data))
model_coef_opt = solution.x

x_test = np.linspace(0,3000,1000)
a_opt = model_coef_opt[:19]
b_opt = model_coef_opt[19:]
plt.scatter(x_data,y_data)
plt.plot(x_test,model(x_test,a_opt,b_opt))

Graphical Fit

Related