[AstroPy] Chi-square problems with lmfit and scipy

[Peter Erwin]

Hi Rudolf,

As I suggested in my second comment, you haven’t actually “solved” your problem
by normalizing the y values; you’ve just transformed one (meaningless) chi^2 value
into a different (meaningless) value.

What I meant by “missing errors” was not “you need to invent some errors and add
them to your data” — which is what your “y = y + y_noise” is doing. I meant that
the data should ideally come with noise estimates (e.g., from whatever pipeline
processed the original data).

If all you have is the y values by themselves (in the form of e.g. flux values, as yours
seem to be), then there’s really nothing to be done. Just be aware that the fit is not
ideal (because all the points are being weighted equally, regardless of their actual
S/N) and the reduced chi^2 value has no useful information about the “goodness”
of the fit.

cheers,

Peter

> On May 20, 2021, at 11:11 AM, Rudolf Baer <rbaer25 at gmail.com> wrote:
> 
> Hi Peter
> thanks for your reply.
> I have solved the problem by normalizing the y values as suggested by Ivan.
> With my first e mail I provided the python code. I attach the actual data file. 
> Your comment on the errors is correct; the file has no errors. I have introduced a gaussian error, 
> y_noise=(0.02*y * np.random.normal(size=len(x)))
> y=y+y_noise
> but it does not affect the results.
> 
> with kind regards
> Rudolf
> 
> 
> On Wed, May 19, 2021 at 10:38 PM Peter Erwin <peter.erwin at gmail.com> wrote:
> 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
> 
> 
> 
> 
> _______________________________________________
> AstroPy mailing list
> AstroPy at python.org
> https://mail.python.org/mailman/listinfo/astropy
> <baseline_BASS0006 0407.csv>_______________________________________________
> 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