Digamma function fitting

[Max Krokever]

Dear All,

I am attempting to fit a version of the Hikami-Larkin-Nagaoka expression to a set of data which involves 2 digamma functions and a natural logarithm:

where

B is the independent variable, the fitting parameters I take to be t_phi, t_e and D, and the rest are constants.

An attempt at a fitting results in a "NaN values detected" type of error (full traceback below). I tried tweaking the starting values of the parameters as well as simplifying the formula. I notice that removing the logarithmic term removes the NaN error but then t_phi and t_e remain stuck at their original values. The whole expression evaluates to approx. 8e-05 and the logarithm evaluates to approx. 6.9, when taking the initial values of the parameters (I had a look at the NaN section in the FAQ).

I would appreciate any suggestions that you may have to deal with the NaN error and/or possibly obtain a fit of the data. I include a link to a file with the data, the code and a full traceback below.

Thank you,

Max

The data:

https://www.dropbox.com/s/9od25ktuvkibmr5/magnetocond.csv?dl=0

The code:

import pandas as pd
import numpy as np
import scipy.constants as sc
import scipy.special as ss
from lmfit import Minimizer, Parameters, report_fit
import matplotlib.pyplot as plt

data=pd.read_csv("magnetocond.csv", header=0)

#take independent and dependent variable in the range of -0.5 to 0.5T
B = np.array(data['Field [T]'])[38:49] #the applied magnetic field in the range of -0.5 to 0.5T
mcond = np.array(data['Conductivity [S/cm]'] - data.iloc[43]['Conductivity [S/cm]'])[38:49] #relative conductivity with respect to the zero magnetic field value


# HLN formula
def HLN(params, B, mcond):
    t_phi = params['t_phi']
    t_e = params['t_e']
    D = params['D']
    model = (sc.e**2/(2*sc.pi**2*sc.hbar))*(ss.digamma(0.5+(sc.hbar/(4*sc.e**t_phi*D*abs(B+0.00001))))-ss.digamma(0.5+(sc.hbar/(4*sc.e*t_e*D*abs(B+0.00001))))+np.log(t_phi/t_e))
    return model - mcond

# Parameters
params = Parameters()
params.add('t_phi', value=1E-9)
params.add('t_e', value=1E-12)
params.add('D', value=1)

# The fitting
result = Minimizer(HLN, params, fcn_args=(B, mcond)).minimize()

# The final result
fit = mcond + result.residual

# Error report
report_fit(result)

# Plot results
try:
    plt.figure(figsize=(10,10))
    plt.plot(B, mcond, 'bo')
    plt.plot(B, fit, 'r')
    plt.xlabel('$\mu_{0}B$ (T)', fontsize=18)
    plt.ylabel('$\sigma - \sigma_{0}$ (S/cm)', fontsize=18)
    plt.show()
except ImportError:
    pass

The full traceback:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-20-14c0825f71c7> in <module>
     28
     29 # The fitting
---> 30 result = Minimizer(HLN, params, fcn_args=(B, mcond)).minimize()
     31
     32 # The final result

~\anaconda3\lib\site-packages\lmfit\minimizer.py in minimize(self, method, params, **kws)
   2361                         val.lower().startswith(user_method)):
   2362                     kwargs['method'] = val
-> 2363         return function(**kwargs)
   2364
   2365

~\anaconda3\lib\site-packages\lmfit\minimizer.py in leastsq(self, params, max_nfev, **kws)
   1698         result.call_kws = lskws
   1699         try:
-> 1700             lsout = scipy_leastsq(self.__residual, variables, **lskws)
   1701         except AbortFitException:
   1702             pass

~\anaconda3\lib\site-packages\scipy\optimize\minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    420         if maxfev == 0:
    421             maxfev = 200*(n + 1)
--> 422         retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol,
    423                                  gtol, maxfev, epsfcn, factor, diag)
    424     else:

~\anaconda3\lib\site-packages\lmfit\minimizer.py in __residual(self, fvars, apply_bounds_transformation)
    602             raise AbortFitException("fit aborted by user.")
    603         else:
--> 604             return _nan_policy(np.asarray(out).ravel(),
    605                                nan_policy=self.nan_policy)
    606

~\anaconda3\lib\site-packages\lmfit\minimizer.py in _nan_policy(arr, nan_policy, handle_inf)
   2445                    'handle this! Please read https://lmfit.github.io/lmfit-py/faq.html#i-get-errors-from-nan-in-my-fit-what-can-i-do '
   2446                    'for more information.')
-> 2447             raise ValueError(msg)
   2448     return arr
   2449

ValueError: NaN values detected in your input data or the output of your objective/model function - fitting algorithms cannot handle this! Please read https://lmfit.github.io/lmfit-py/faq.html#i-get-errors-from-nan-in-my-fit-what-can-i-do for more information.

[Matt Newville]

On Wed, Jul 14, 2021 at 1:01 PM Max Krokever <max.kr...@gmail.com> wrote:

Dear All,

I am attempting to fit a version of the Hikami-Larkin-Nagaoka expression to a set of data which involves 2 digamma functions and a natural logarithm:

where

B is the independent variable, the fitting parameters I take to be t_phi, t_e and D, and the rest are constants.

An attempt at a fitting results in a "NaN values detected" type of error (full traceback below). I tried tweaking the starting values of the parameters as well as simplifying the formula. I notice that removing the logarithmic term removes the NaN error but then t_phi and t_e remain stuck at their original values. The whole expression evaluates to approx. 8e-05 and the logarithm evaluates to approx. 6.9, when taking the initial values of the parameters (I had a look at the NaN section in the FAQ).

If you read the FAQ section, then you did read under the section entitled "Common sources of NaN" that "taking sqrt(x) or log(x) where x is negative" is listed as the first, most common source of NaN.  You even add that removing the logarithmic term removes the NaN.

So, what would happen in your code if `t_phi` or `t_e` were negative?  Yes, that's correct, your function will give NaN.  As you read in the FAQ, a NaN that comes from your calculation in your model or objective function cannot be fixed by the fitting code or by "us".  It has to be fixed by you.  Are you doing anything to ensure that the value of `t_phi / t_e` cannot ever be negative? I believe you are not.

I would appreciate any suggestions that you may have to deal with the NaN error and/or possibly obtain a fit of the data. I include a link to a file with the data, the code and a full traceback below.

All of the advice I am able to give is in the FAQ.  If you are taking the log or square root of some value, you really have to ensure that value is positive.  There is no way around this and no way to make the advice be more encouraging, enabling, or user-friendly than saying "read the FAQ". 

[Max Krokever]

I would think that the appearance of NaNs can also depend on how the initial parameter values are varied (this is what I tried to imply through the starting values not causing problems). Without an understanding of how this is done here I will resort to making sure there is no chance of a NaN appearing for any values across the expression, as you sufficiently emphasize. Thank you for your time.

[Marc William Pound]

Hello Max,

You can set the range of allowed values for Parameters with min and max keywords to params.add().  The default for these are -inf and inf, so it is usually a good idea to limit the search range.   This may help with your NaN problem if it is coming from the log term and it is valid for your problem to limit t_e and t_phi to positive values.

best,
Marc