CODEWITHSUNDEEP
optimizationscipycurve-fitting
Pete Lavelle
Pete Lavelle

Reputation: 103

Curve Fitting in Python using scipy.optimize curve_fit

As a relative beginner in Python, i'm struggling to understand (and therefore use) the "curve_fit" from scipy.optimize

I've tried following answers to previous questions: python numpy/scipy curve fitting and exponential curve fitting with python but unfortunately had no luck getting it to work.

This is an engineering problem (dealing with measured data from a test rig), therefore I know that the formula is in the format Y=(X/A)+(X/B)^(1/C) where A,B,C are the constants that need to be found

[Current Code]

import numpy as np
from scipy.optimize import curve_fit

valY_list = [yyy1,yyy2,yyy3,yyy4]
valX_list = [xxx1,xxx2,xxx3,xxx4]

val_Y = np.array(valY_list)
val_X = np.array(valX_list)

def fit_func(val_X,A,B,C):
    return (val_X/A)+((val_X/B)^(1/C))

params = curve_fit(fit_func, val_X, val_Y)

[A,B,C] = params[0]

NB1: in reality valY_list & valX_list are >500 entries long (stored as floats).

NB2: I also know that the values of A, B, C should be within a certain range of values, so I want to constrain the solution when it performs the optimisation.

0.005 < A < 0.5

0.0 < B < 5000

0.0001 < C < 0.1

I realise that my code is probably quite rudimentary, and likely has lots of things missing (or glaringly obvious mistakes for an experienced coder!) so my apologies. Any help would be appreciated!

Upvotes: 2

Views: 4211

Answers (1)

M Newville
M Newville

Reputation: 7862

You might find lmfit (https://lmfit.github.io/lmfit-py) useful for this and other curve-fitting problems. It provides a higher-level interface to curve-fitting than curve_fit and has many convenient and advanced options for model building and working with parameters and fit statistics.

With lmfit, I would suggest the following approach:

import numpy as np
from scipy.optimize import curve_fit
from lmfit import Model

valY_list = [yyy1,yyy2,yyy3,yyy4]
valX_list = [xxx1,xxx2,xxx3,xxx4]

val_Y = np.array(valY_list)
val_X = np.array(valX_list)

def fit_func(x, a, b, c):
    return (a*x)+(b*x)**c

mymodel = Model(fit_func)

params = mymodel.make_params(a=10, b=1, c=100.0)

params['a'].min = 2.0
params['a'].max = 200.

params['b'].min = 0.002
params['b'].max = 1.e8

params['c'].min = 10.0
params['c'].max = 10000.0

result = mymodel.fit(val_Y, params, x=val_X)

print(result.fit_report())

for par_name, param in result.params.items():
    print(par_name, param.value, param.stderr)

Note that I changed your model function slightly to use 1/PARAM, and adjusted the bounds accordingly (I think!). The printed fit report will include fit statistics and the best fit values and standard errors for each of the variables. Also, note that with lmfit that Parameters are named according to the argument names of your fitting function, and that min/max bounds go with the Parameter object. Although the example above sets bounds, you could also fix any of the Parameters, with (for example):

params['c'].vary = False

And, if you really do want your A, B, and C, you can make parameters for these as constrained expressions:

params.add('A', expr='1/a')
params.add('B', expr='1/b')
params.add('C', expr='1/c')

These won't vary in the fit, but their values and standard errors will be reported. Essentially any valid python expression, using other parameter names and basic math functions can be used.

There are many more features and decent documentation and examples, but that should get you started.

Upvotes: 1

Related Questions