[AstroPy] Chi-square problems with lmfit and scipy

[Peter Erwin]

Hi Rudolf,

My apologies -- I had failed to see that you *had* included the code you used in
your original post.
So, please ignore my previous email.


That said (and having taken a closer look at lmfit), it's clear that the "chi^2"
value returned by lmfit in your example is almost certainly what you *should* get,
given that you have not provided any error values to go with your data. So it's
almost certainly not some effect of "numerical accuracy limits".

The Model.fit method of lmfit has an optional input called "weights". For
standard chi^2 minimization, these should have values of 1/sigma, where sigma
are the error values corresponding to the data values, *in the same units*.
If you *don't* supply a "weights" vector, then lmfit effectively treats all
data points have sigma = 1, which is of course not true for your data.

What lmfit minimizes is the sum of (weights*(data - model))^2. If weights = 1/sigma,
then this is the sum of ((data - model)/sigma)^2, which is the standard chi^2
formula for data with errors, and the "reduced" version of this will (assuming the 
errors are accurate and quasi-Gaussian) approximate a chi^2 distribution
and can be used (with suitable caution) as a crude "goodness of fit" indicator.

But if weights=None [the default if you don't supply a weights keyword], then 
internally lmfit minimizes the sum of (data - model)^2.
The "reduced chi^2" version of this is simply [sum of (data - model)^2]/(N_data - N_var),
which for large N_data is approximately the mean value of (data - model)^2.

You don't provide any examples of data values, I'm guessing they're somewhere
in the neighborhood of 10^11 or 10^-12. This is because that's the square root of your
"reduced chi^2" value, and is something like the RMS value of (data - model).


So it's important to understand that without providing the correct error values
as input to the Model.fit method (via the "weights" keyword), the "chi^2" and
"reduced chi^2" values returned by lmfit are *meaningless*. Rescaling your data
so it has maximum values of 1 will give you different but *equally meaningless*
values of chi^2 and reduced chi^2.

cheers,

Peter


> On May 19, 2021, at 5:43 PM, Rudolf Baer <rbaer25 at gmail.com> wrote:
> 
> I am fitting the continuum of a NIR spectrum with a blackbody plus power law using lmfit and scipy. The plot of the fit looks good
> and the resulting parameters agree reasonably with previous results.
> T:      1109.29505 +/- 3.08340774 (0.28%) (init = 1100)
> const:  4.7240e-21 +/- 6.3014e-23 (1.33%) (init = 2e-21)
> A:      2.4188e-11 +/- 6.3181e-14 (0.26%) (init = 2e-11)
> p:     -1.06704890 +/- 0.01701493 (1.59%) (init = -1)
> 
> However the chi-square and reduced chi-square values
> (chi-square = 1.1268e-20 reduced chi-square = 2.1039e-24)
> are clearly wrong. The problem is apparently known as the scipy guide says Note that the calculation of chi-square and
> reduced chi-square assume that the returned residual function is scaled properly to the uncertainties in the data. 
> For these statistics to be meaningful, the person writing the function to be minimized must scale them properly.
> 
>  I do not know how to do this. Any advice will be appreciated.
> 
> With kind regards
> Rudolf Baer
> Zurich
> 
>  
> 
> The input data for the code are x: lambda in AA, y: flux density in erg/cm^{2}/s/{\AA)
> Code:
> 
> x=x_AA/10000. #<---------conversion AA to um
> 
> y=y*1e4
> 
> y=y*x          
> #=========================================================================================
> 
> def bb(x, T, const):
> 
>     from scipy.constants import h,k,c
> 
>     x = 1e-6 * x # convert to metres from um
> 
>     return const*2*h*c**2 / (x**5 * (np.exp(h*c / (x*k*T)) - 1)) #J/s/m2/m
> 
>  
> def powerlaw(x,A,p):
> 
>     return A*x**p
> 
>  
> mod= Model(bb) + Model(powerlaw)
> 
> pars  = mod.make_params(T=1100,const=2*1e-21,A=2*1e-11,p=-1.0)
> 
>  result = mod.fit(y,pars,x=x)
> _______________________________________________
> AstroPy mailing list
> AstroPy at python.org
> https://mail.python.org/mailman/listinfo/astropy

=============================================================
Peter Erwin                   Max-Planck-Insitute for Extraterrestrial 
erwin at mpe.mpg.de              Physics, Giessenbachstrasse
tel. +49 (0)176 2481 7713     85748 Garching, Germany
fax  +49 (0)89 30000 3495     https://www.mpe.mpg.de/~erwin