How to fit log(a-x) type functions with scipy. curve_fit?

user_8932869 0
Viewed 752

I am trying to fit a function which looks like log(y)=a*log(b-x)+c, where a, b and c are the parameters that need to be fitted. The relevant bit of code is

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

def logfunc(T, a, b, c):
    v=(a*np.log(b-T))+c
    return v

popt, pcov=curve_fit(logfunc, T, np.log(Energy), check_finite=False, bounds=([0.1, 1.8, 0.1], [1.0, 2.6, 1.0]))

plt.plot(T, logfunc(T, *popt))
plt.show

Where T and Energy is some data that was generated (I use it to plot other things so the data should be fine). T is between 0.3 and 3.2. I am pretty sure that the problem is the fact that there is a point where b=T because I keep getting the error ValueError: Residuals are not finite in the initial point. but I am not sure how to solve this.

curve-fitting scipy python-3.x numpy
2 Answers

You may find the lmfit package (http://lmfit.github.io/lmfit-py/) useful for this sort of problem. This provides a higher-level approach to curve fitting problems and a better abstraction of Parameters and Models than scipy.optimize package or curve_fit() function.

For the problem here, two important features of lmfit are

  1. the ability to set bounds on variables. curve_fit() can do this as well, but only by working with ordered lists of min/max bounds. With lmfit, the bounds belong to Parameter objects.
  2. having a way to explicitly set a policy for handling NaN values, which could definitely cause problems for your fit.

With lmfit, your script would be written approximately as

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model

def logfunc(T, a, b, c):
    return (a*np.log(b-T))+c

log_model = Model(logfunc, nan_policy='raise')  # raise error on NaNs
params = log_model.make_params(a=0.5, b=2.0, c=0.5) # initial values
params['b'].min = 1.8  # set min/max values
params['b'].max = 2.6 
params['c'].min = 0.1  # and so forth 

result = log_model.fit(np.log(Energy), params, T=T)

print(result.fit_report())

plt.plot(T, Energy, 'bo', label='data')
plt.plot(T, np.exp(result.best_fit), 'r--', label='fit')
plt.legend()
plt.xlabel('T')
plt.ylabel('Energy')
plt.gca().set_yscale('log', basey=10)
plt.show()

This is slightly more verbose than your starting script because it gives a labeled plot and because using Parameter objects instead of scalars gives more flexibility and clarity.

For your fit, you might consider setting the nan_policy to 'omit', which will omit NaNs as they occur -- never a great idea, but sometimes helpful to get you started on finding where log(b-T) is valid. You could also alter your model function to do something like

def logfunc(T, a, b, c):
    arg = b - T
    arg[np.where(arg < 1.e-16)] = 1.e-16
    return a*np.log(arg) + c

To explicitly prevent one obvious cause of NaNs.

user_51797480

Residuals are not finite in the initial point

means the initial point is bad, where some logarithms are infinite or undefined. You need a better initial point.

By the nature of the model, b has to be greater than any of the points in T. The bounds on b that you have at present do not guarantee that. Tighten them up.

When you do not provide p0 parameter, SciPy will take a guess within the provided bounds. So if the bounds guarantee finiteness, the error will not occur. Still, it is generally better to prescribe p0 yourself, because you have better a priori understanding of the problem than SciPy does.

A working example with adjusted bounds:

popt, pcov=curve_fit(logfunc, np.linspace(0.3, 3.2, 6), [8, 7, 6, 5, 4, 3], bounds=([0.1, 3.2, 0.1], [1.0, 3.6, 1.0]))
Admin1
Related