about exponential fitting

[Tetuko Kurniawan]

Hi all and Matt Newville,

I'm beginner, and I'm doing a very basic thing about exponential fitting.

However, i have some difficulties that i couldn't find the answer, although I try many suggestions from previous post or in stackoverflow.

Here is the model:

def mod_exp(x, amplitude, decay, x0, y0):
       return amplitude*np.exp( -(x-x0)/decay)+y0

 mod = lt.Model(mod_exp)
 mod.set_param_hint('amplitude', value=0.8, vary=True, min =0.0)
 mod.set_param_hint('decay', value=2, vary=True, min=0.1)
 mod.set_param_hint('x0', value=4, vary=True, min=0.0)
 mod.set_param_hint('y0', value=0, vary=True)
 pars = mod.make_params()
 out = mod.fit(eps, pars, x=x_eps)

The example data is attached. X is the independent variable.

If you do the data with above example, you will get the result quite ok..but:

1. As you might notice, there's a possibility the -(x-x0)/decay, becomes positive and very big. First if 'decay' too small and if the (x-x0) coincidentally becomes <0.

if let say change value of 'decay' = 1, you will get the nan_policy warning.

if i modified

mod.set_param_hint('x0', value=4, vary=True, min=np.min(X))

It prevents negative value of x-x0.

But the fit result of x0 stuck on np.min(X)... I think it is not the best fitting, although it works better.

example  of fit report:

Qc=4.0,Qd=1.6
[[Model]]
    Model(mod_exp)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 26
    # data points      = 30
    # variables        = 4
    chi-square         = 0.01068116
    reduced chi-square = 4.1081e-04
    Akaike info crit   = -230.214133
    Bayesian info crit = -224.609343
##  Warning: uncertainties could not be estimated:
    x0:         at initial value
    x0:         at boundary
[[Variables]]
    amplitude:  0.53622335 (init = 0.8)
    decay:      0.38896822 (init = 1)
    x0:         4.16521739 (init = 4.165217)
    y0:         0.29091852 (init = 0)

2. using python 'min' as suggested from other example:

amplitude*np.exp(min(700, -(x-x0)/decay))+y0

i got this error:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all().

3. I tried using:

mod.set_param_hint('decay', value=1, vary=True, min=np.min(x_eps)/700)

I think x0 should be counted on the minimum value, but x0 is another parameter, I don't know how to do it.

I'm not sure about this. Any suggestion?

I'm sorry if I this is trivial/beginner question, but I want to understand about this.

Thank you very much!!

Best regards,

Teko

[Matt Newville]

See responses below:

The fit report is saying that "x0" is both at the initial value and the boundary.   You had

       mod.set_param_hint('x0', value=4, vary=True, min=np.min(X))

and the minimum value of the X array you pass in is 4.165217.   

First, do not give initial values that exceed the bounds.  That just makes no sense.   Bounds mean that the Parameter cannot ever have values outside of those limits.  It seems to me that about half of the "Help, I got a bad fit" is the result of poor use of bounds and setting them either arbitrarily too tight or automatically set based on input data.  Use bounds sparingly and carefully.

Second, don't use `set_param_hint()` that way.  I am continuously dumbfounded by this usage.  Parameter hints belong **TO THE MODEL**.  Is it a property of the model that the initial value for `x0` is 4?  No, it is not.  That is a property of the data.   Do not use `set_param_hint()` this way.  

2. using python 'min' as suggested from other example:

amplitude*np.exp(min(700, -(x-x0)/decay))+y0

i got this error:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all().

Yes, that is how Python's `min()` works.  You might want to read about `numpy.where`.

3. I tried using:

mod.set_param_hint('decay', value=1, vary=True, min=np.min(x_eps)/700)

I think x0 should be counted on the minimum value, but x0 is another parameter, I don't know how to do it.

Again, setting parameter hints in this way is an exceptionally bad idea.  Automatically setting bounds based on other values is a very easy way to get an error, as you got above.

I'm not sure about this. Any suggestion?

I'm sorry if I this is trivial/beginner question, but I want to understand about this.

Thank you very much!!

Other advice:   If there is an exponential decay to your data and model, you might find it helpful fit the log of your data to a model that takes the log of whatever your model currently returns.

Best regards,

Teko

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/e2f848e2-f70e-41fe-a5dc-efb697a672f7n%40googlegroups.com.

--Matt 

[Tetuko Kurniawan]

Hi Matt,

Thank you for your reply, and also I would to say thank you for making this library available for everyone.

Recently I've been using this more often than Scipy optimize, but I have to learn more on this.

One key point you mention is about using the parameter objects rather than 'set_param_hint()'.

I'll keep in mind about that while reading again your documentation.

I was using 'set_param_hint()' because it seems shorter, and (i thought) it works the same.

The value= 4 was stupid, I  miss that.

Yes, the log idea should work better.

I'll keep working on this.

Thank you.

Teko

[Tetuko Kurniawan]

Hi Matt,

I hope you have time for this.

Actually I'm going to use the log, but I'm still curious about 'set_param_hint'.

So i modified the script, it is not perfect, but works better than before:

###

def mod_exp(x, amplitude, decay, x0, y0):

    if np.all(-(x-x0)/decay) >= 500:
        arg = 500
    else:
        arg = -(x-x0)/decay
    return amplitude*np.exp(arg)+y0

mod = lt.Model(mod_exp)
params = mod.make_params()
params['amplitude'].min = 0.0
params['amplitude'].vary = True
params['amplitude'].value = 0.8
params['decay'].min = 0.0
params['decay'].vary = True
params['decay'].value = 5
params['x0'].min = 0.0
params['x0'].vary = True
params['x0'].value = 5
params['y0'].vary = True
params['y0'].value = 0

out = mod.fit(Y, params=params, x=X)  # Using the same data as before but i provide it above

for pname, par in params.items():
        print(pname, par)
print(out.fit_report(min_correl=0.25))

###

Output:

amplitude <Parameter 'amplitude', value=0.8, bounds=[0.0:inf]>
decay <Parameter 'decay', value=5, bounds=[0.0:inf]>
x0 <Parameter 'x0', value=5, bounds=[0.0:inf]>
y0 <Parameter 'y0', value=0, bounds=[-inf:inf]>

[[Model]]
    Model(mod_exp)
[[Fit Statistics]]
    # fitting method   = leastsq

    # function evals   = 71

    # data points      = 30
    # variables        = 4
    chi-square         = 0.01068116
    reduced chi-square = 4.1081e-04

    Akaike info crit   = -230.214130
    Bayesian info crit = -224.609341
[[Variables]]
    amplitude:  0.54158662 +/- 691538.645 (127687541.84%) (init = 0.8)
    decay:      0.38896214 +/- 0.02346960 (6.03%) (init = 5)
    x0:         4.16133832 +/- 499461.115 (12002415.47%) (init = 5)
    y0:         0.29092002 +/- 0.00769647 (2.65%) (init = 0)
[[Correlations]] (unreported correlations are < 0.250)
    C(amplitude, x0) = -1.000
    C(decay, y0)     = -0.747

as above create parameters object that has name params, which consist of 4 parameter.

if i refer back to previous code:

    mod.set_param_hint('amplitude', value=0.8, vary=True, min=0.0)

    mod.set_param_hint('decay', value=5, vary=True, min=0.0)
    mod.set_param_hint('x0', value=5, vary=True, min=0.0)

    mod.set_param_hint('y0', value=0, vary=True)

    params = mod.make_params()

It is also creating the same parameters object.

Which part should I take a look to see the difference, as you said in the previous email?

Or both ways are still incorrect?

Thank you!

Best,

Teko