scipy.optimize error power function

Question

I have some trouble to fit a set of value with the given function : f(x)= const*(1-(x/a)**b)**c

I am using python 3.6.3 with the following code :

import numpy as np
import scipy.optimize as opt
from scipy.optimize import curve_fit

x=[0.,0.4,0.8,1.6,2.,2.4]
y=[0.09882902,0.07298427,0.05111438,0.01679405,0.00517385,0.00065633]

def func(x,a,b,c):
    return y[0] * ( 1 - (x/a)**b )**c

x0=np.array([2.0,0.9,1.5])
opt.curve_fit(func,x,y,p0=x0)

I have the following error message:

RuntimeWarning: invalid value encountered in power
  return y[0] * ( 1 - (x/a)**b )**c
///: OptimizeWarning: Covariance of the parameters could not be estimated
  category=OptimizeWarning)

The problem seems to appear because of the last exponent, because the following function will work fine :

def func(x,a,b,c):
    return y[0] * ( 1 - (x/a)**b )*c

sascha · Accepted Answer · 2017-11-22T16:20:07.220

2

Just do the classic debugging steps: print your components to check if math-ops are well-defined.

Modify your function to:

def func(x,a,b,c):
    print(x/a)
    print((x/a)**b)
    print((1-(x/a)**b))

    result = y[0] * ( 1 - (x/a)**b )**c
    print(result)

    if not np.isfinite(result):
        assert False
    else:
        return result

and you will see:

[ 0.   0.2  0.4  0.8  1.   1.2]
[ 0.          0.23492379  0.43838329  0.81805215  1.          1.17831965]
[ 1.          0.76507621  0.56161671  0.18194785  0.         -0.17831965]
...-py:13: RuntimeWarning: invalid value encountered in power
  result = y[0] * ( 1 - (x/a)**b )**c
[ 0.09882902  0.06613655  0.04159532  0.00767017  0.                 nan]

where the exponentiation does not work for that negative value and a nan is introduced (which will probably introduce many more in later steps).

Remark: initial-point should be all-ones by default (easy to print out too).

Fixing this depends on what you are actually trying to do (change of model; using bounds, ...).

edited Nov 22 '17 at 16:20

answered Nov 22 '17 at 16:14

sascha

32,238
6
68
110

Thank you. Actually, this function worked fine with gnuplot fitting, and I am supposed to use it, though it becomes infinitly small for x -> a If I am not mistaken, the "NaN" seems to be a number below 1e-13. Is there a way to overcome this limitation ? – ant.kr Nov 22 '17 at 16:53
The NaN is **Not a Number** :-) Again: this is all depending on many model-decisions you have access to. If gnuplot worked, than there are differences in the assumptions or the data as the math-problem won't go away (for the general case). – sascha Nov 22 '17 at 16:55
Maybe [this is relevant to you](https://stackoverflow.com/a/45384691/2320035) (and maybe a difference to gnuplot). – sascha Nov 22 '17 at 17:00
Well, I know what nan means. If I run the curve fitting, starting from [a=2.4,b=0.9,c=1.5], the value of func(x) becomes "not a number" just after a step where it was equal to 1e-13. So maybe the number of digits becomes too important? – ant.kr Nov 22 '17 at 17:05
Becoming a nan there has nothing to do with a threshold. – sascha Nov 22 '17 at 17:07
Well, my mistake. I will try with some bounds – ant.kr Nov 22 '17 at 17:25

score 0 · Answer 2 · answered Nov 23 '17 at 16:16

I found the solution, thanks Sascha. The function is indeed undetermined when x > a. To overcome this problem, I constrained the value of a to be always superior to the max value of x :

fit=opt.curve_fit(func,x,y,p0=x0,bounds=([x[-1]+0.001,0.0,1.00],[5.0,1.0,2.0]))

However, in gnuplot, the fit was possibly finding x > a. I dont know why. Maybe it is taking the real part of func(x>a), but I dont know if it is really straightforward to do so.

scipy.optimize error power function

2 Answers2