Scaling error and Information cirteria in lmfit

Question

I am using lmfit for solving a non-linear optimization problem. It works fine to the point where I am trying to implement a measurement error as the standard deviation of the dependent variable y (sigma_y). My problem is, that I cannot interpret the Information criteria properly. When implementing the return (model - y)/(sigma_y) they just raise from really low to very high values.

i.e. [left: return (model - y) -weighting-> right: return (model - y)/(sigma_y)]:

• chi-square 0.00159805 -> 47.3184972
• reduced chi-square 1.7756e-04 -> 5.25761080 expectation value is 1 || SO discussion
• Akaike info crit -93.2055413 -> 20.0490661 the more negative, the better
• Bayesian info crit -92.4097507 -> 20.8448566 the more negative, the better

My guess is, that this is somehow connected to bad usage of lmfit (wrong calculation of Information criteria, bad error scaling) or to a general lack of understanding these criteria (to me reduced chi-square of 5.258 (under-estimated) or 1.776e-4 (over-estimated) sounds like a really poor fit, but the plot of residuals etc. looks okay for me...)

Here is my example code that reproduces the problem:

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

y = np.array([1.42774324, 1.45919099, 1.58891696, 1.78432768, 1.97403125,
       2.17091161, 2.02065766, 1.83449279, 1.64412613, 1.47265405,
       1.4507    ])
sigma_y = np.array([0.00312633, 0.00391546, 0.00873894, 0.01252675, 0.00639111,
       0.00452488, 0.0050566 , 0.01127185, 0.01081748, 0.00227587,
       0.00190191])
x = np.array([0.02372331, 0.07251188, 0.50685822, 1.30761963, 2.10535442,
       2.90597497, 2.30733552, 1.50906939, 0.7098041 , 0.09580686,
       0.04777082])
offset = np.array([1.18151707, 1.1815602 , 1.18202847, 1.18187962, 1.18047493,
       1.17903493, 1.17962602, 1.18141625, 1.18216907, 1.18222051,
       1.18238949])
parameter = lmfit.Parameters()
parameter.add('par_1', value = 1e-5)
parameter.add('par_2', value = 0.18)

def U_reg(parameter, x, offset):
    par_1 = parameter['par_1']
    par_2 = parameter['par_2']
    model = offset + 0.03043211217 * np.arcsinh(x / (2 * par_1)) + par_2 * x
    return (model - y)/(sigma_y)

mini = lmfit.Minimizer(U_reg, parameter, fcn_args=(x, offset), nan_policy='omit', scale_covar=False)
regr1 = mini.minimize(method='least_sq') #Levenberg-Marquardt
print("Levenberg-Marquardt:")
print(lmfit.fit_report(regr1, min_correl=0.9))

def U_plot(parameter, x, offset):
    par_1 = parameter['par_1']
    par_2 = parameter['par_2']
    model = offset + 0.03043211217 * np.arcsinh(x / (2 * par_1)) + par_2 * x
    return model - y

plt.title("Model vs Data")
plt.plot(x, y, 'b')
plt.plot(x, U_plot(regr1.params, x, offset) + y, 'r--', label='Levenberg-Marquardt')
plt.ylabel("y")
plt.xlabel("x")
plt.legend(loc='best')
plt.show()

I know that it might be more convient to implement weighting via lmfit.Model, due to multiple regressors I want to keep the lmfit.Minimizer method. My python version is 3.8.5, lmfit is 1.0.2

Answer 1

Well, in order for the magnitude of chi-square to be meaningful (for example, that it be around (N_data - N_varys), the scale of the uncertainties has to be correct -- giving the 1-sigma standard deviation for each observation.

It's not really possible for lmfit to detect whether the sigma you use has the right scale.