[AstroPy] Chi-square problems with lmfit and scipy

Fri May 21 04:22:13 EDT 2021

Hi  Peter
thank you very much for your help.
Cheers
Rudolf

On Thu, May 20, 2021 at 2:37 PM Peter Erwin <erwin at mpe.mpg.de> wrote:

> 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
>
>
>
>
> _______________________________________________
> AstroPy mailing list
> AstroPy at python.org
> https://mail.python.org/mailman/listinfo/astropy
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mail.python.org/pipermail/astropy/attachments/20210521/b30e7ad1/attachment-0001.html>