[AstroPy] Chi-square problems with lmfit and scipy

[Peter Erwin]

Rudolf,

It’s not clear to me what you are doing. You say you are using “lmfit and scipy”.
By “lmfit” do you mean the lmfit package — https://lmfit.github.io/lmfit-py/index.html ?

This has nothing to do with scipy (though it assumes numpy exists), so maybe you
are using something else instead?

In any case, I’d guess the chi^2 values are “wrong” because you don’t have good error
estimates (“uncertainties”) for the data values.

Providing examples of the actual Python commands you are using (and maybe a few
sample data values from your NIR spectrum) would be helpful.

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