scipy.odr - Goodness of fit and parameter estimation for explicit orthogonal distance regression
Hi list, Currently, I am trying to fit a quadratic curve to a data set which has much larger errors in the x than in the y direction. My errors are assumed to be normally distributed and I want to estimate the confidence interval of the fitted parameters. I have fitted the data two different ways. 1) I neglect the x errors and fit the quadratic by minimizing the weighted residuals (y -f) / sig_y via scipy.optimize.leastsq and 2) I use scipy.odr to fit the parameters. Both result similar fitted parameters. Now I am stuck with estimating the confidence intervals on these errors and I have a couple of questions. For method 1) the reduced chi squared is bad (much larger then 1), because when neglecting the x errors, none of the points actually lie within a couple of standard deviations on the line. However when including the x-errors they all fall nicely onto the line. Is there a way for me to include the x errors into my minimization and goodness of fit estimation via the reduced chi-squared? I came across a remark in the CERN minuit documentation [1], that seemed to approximate the function by a line over the point and then use this to convert x to y errors. I also read something like this in numerical recipes [2], but not for general least squares (only for the case of linear data). Does anyone have any pointers into that direction? My second question is related to method 2). Is there a way of accessing the goodness of fit for ODR, similar to calculating the reduced chi-squared for a fit that only has errors in the y? Another question along this line concerns the scipy.odr.ODR.Output.sd_beta attribute. In the docstring it says it is the standard error of the parameter, does that mean a 1 standard deviation confidence interval? And how exactly are they calculated. Tried to look at the source and in the odrpack guide, but unfortunately couldn't figure that out. I somehow have the feeling, that my problem should be quite standard (goodness of fit with both x and y errors), but so far I could not find a good explanation. Any pointers to textbooks or resources on the web would be greatly appreciated. Best regards, Markus [1] http://wwwasd.web.cern.ch/wwwasd/cgi-bin/listpawfaqs.pl/27 [2] Numerical Recipes in C, Second Edition, Chapter 15.3 "Straight-Line Data with Errors in Both Coordinate"
Perhaps the module kmpfit from the kapetyn package would be helpful to you. http://www.astro.rug.nl/software/kapteyn/kmpfit.html Don't miss the tutorial http://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html which has a section "Fitting data when both variables have uncertainties" with subsection "Effective variance method for various models: Model parabola"
Don't miss the tutorial http://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html which has a section "Fitting data when both variables have uncertainties"
Thanks a lot, the tutorial and references therein are very useful. Regards, Markus
On Mon, May 14, 2012 at 4:25 AM, Markus Baden <markus.baden@gmail.com> wrote:
Hi list,
Currently, I am trying to fit a quadratic curve to a data set which has much larger errors in the x than in the y direction. My errors are assumed to be normally distributed and I want to estimate the confidence interval of the fitted parameters. I have fitted the data two different ways. 1) I neglect the x errors and fit the quadratic by minimizing the weighted residuals (y -f) / sig_y via scipy.optimize.leastsq and 2) I use scipy.odr to fit the parameters. Both result similar fitted parameters.
Now I am stuck with estimating the confidence intervals on these errors and I have a couple of questions.
scipy.odr provides an estimate of the covariance matrix and standard deviations of the parameter estimates. Getting the confidence interval for a parameter is just a matter of scaling up the standard deviations by the appropriate t-distribution value with nobs-nparams degrees of freedom. A paper by the ODRPACK implementors gives the formula explicitly on page 6: http://www.mechanicalkern.com/static/odr_vcv.pdf It also has more information on how the covariance matrix is calculated.
My second question is related to method 2). Is there a way of accessing the goodness of fit for ODR, similar to calculating the reduced chi-squared for a fit that only has errors in the y?
Output.res_var is the reduced Chi-square. But you can compute it from scratch using the raw residuals. Output.eps are the differences in Y and Output.delta are the differences in X. Square each and multiply by their weights (a.k.a. divide by the data variances), then add up all of them. Divide by (nobs - nparams).
Another question along this line concerns the scipy.odr.ODR.Output.sd_beta attribute. In the docstring it says it is the standard error of the parameter, does that mean a 1 standard deviation confidence interval?
Yes.
And how exactly are they calculated. Tried to look at the source and in the odrpack guide, but unfortunately couldn't figure that out.
As given in the odr_vcv.pdf paper. Some slightly sparser details are also in the ODRPACK Guide in section "4.B. Covariance Matrix". -- Robert Kern
scipy.odr provides an estimate of the covariance matrix and standard deviations of the parameter estimates. Getting the confidence interval for a parameter is just a matter of scaling up the standard deviations by the appropriate t-distribution value with nobs-nparams degrees of freedom. A paper by the ODRPACK implementors gives the formula explicitly on page 6:
http://www.mechanicalkern.com/static/odr_vcv.pdf
It also has more information on how the covariance matrix is calculated.
Thanks for the pointer to the paper. That is exactly what I was looking for. Thanks also for the other clarification (and making scipy.odr!) Regards, Markus
participants (3)
-
John Gleeson -
Markus Baden -
Robert Kern