MathGroup Archive 2013

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Error in slope and intercept

  • To: mathgroup at
  • Subject: [mg131781] Re: Error in slope and intercept
  • From: Ray Koopman <koopman at>
  • Date: Thu, 3 Oct 2013 00:36:00 -0400 (EDT)
  • Delivered-to:
  • Delivered-to:
  • Delivered-to:
  • Delivered-to:

On Wed, Oct 02, 2013 @ 02:53 AM, Gebbo <nicolasgebbo at> wrote:
> I was informing myself about linear regression with
> error in both observables and found this algorithm:
> (Weighted Orthogonal Regression by Alternating Least Squares)
> worals[x_, y_, sx_, sy_] := Block[
> {a,b,f,z, u = 1/sx, v = 1/sy, w = (sy/sx)^2},
> {a,b} = (y*v).PseudoInverse@{v,x*v}; f = #.#&[(a+b*x-y)v];
> While[f > (z = (x*w + (y-a)b)/(b^2 + w);
>            {a,b} = (y*v).PseudoInverse@{v,z*v};
>            f = #.#&@Join[(z-x)u,(a+b*z-y)v])];
> {f,{a,b}}]
> which gives me {chisquare, {intercept, slope}} as output.
> This works prefectly fine but I'd like to get the error on
> slope and intercept. I don't understand enough about Mathematica
> or the theory behind this algorithm, so i would pref the solution
> or hinds our i can write it myself.
> Thanks Gebbo

The model is {x = z + d, y = a + b*z + e}, where d & e are random errors, 
and a, b, & z are unknowns for which values are to be found that minimize 
the weighted sum of squares f = #.# & @ Join[(z-x)/sx,(a+b*z-y)/sy]. 
We start with z = x, then minimize f alternately with respect to either 
{a,b} or z with the other held constant, until f no longer changes.

To estimate the error in {a,b}, I usually jackknife the solution.

n = Length@x; {f,ab} = worals[x,y,sx,sy]; 
{jab,jc} = {ab + (n-1)(ab-Mean@#), (n-1)^2/n Covariance@#}& @ 

jab is the jackknifed estimate of {a,b}, and jc is the jackknifed estimate 
of its covariance matrix; Sqrt@Diagonal@jc gives the estimates of the two 
standard errors.

See also

  • Prev by Date: Re: Tick Values on Gauges
  • Next by Date: Re: Tick Values on Gauges
  • Previous by thread: Error in slope and intercept
  • Next by thread: Error on slope and intercept