[SciPy-User] scipy.odr - Goodness of fit and parameter estimation for explicit orthogonal distance regression

[Robert Kern]

On Mon, May 14, 2012 at 4:25 AM, Markus Baden <markus.baden at 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