[AstroPy] Chi-square problems with lmfit and scipy

[Rudolf Baer]

Hi Ivan
brilliant ! it seems to work.
Thank you very much
Rudolf


T:      1109.28599 +/- 3.08367329 (0.28%) (init = 1100)
const:  1.0685e-10 +/- 1.4258e-12 (1.33%) (init = 1)
A:      0.54709366 +/- 0.00142887 (0.26%) (init = 1)
p:     -1.06700052 +/- 0.01701320 (1.59%) (init = -1)


# fitting method   = leastsq
    # function evals   = 36
    # data points      = 5360
    # variables        = 4
    chi-square         = 5.76471199
    reduced chi-square = 0.00107631


On Wed, May 19, 2021 at 10:19 PM Ivan Valtchanov <ivvv68 at gmail.com> wrote:

> Hi Rudolf,
>
> I would recommend to normalise the y-values within [0,1], I suspect
> the results you see are due to numerical accuracy limits.
>
> Cheers,
> Ivan V
>
>
> On Wed, 19 May 2021 at 17:43, 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
> _______________________________________________
> 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/20210520/0b8f33c5/attachment.html>