Re: Error in slope and intercept
- To: mathgroup at smc.vnet.net
- Subject: [mg131781] Re: Error in slope and intercept
- From: Ray Koopman <koopman at sfu.ca>
- Date: Thu, 3 Oct 2013 00:36:00 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- Delivered-to: l-mathgroup@wolfram.com
- Delivered-to: mathgroup-outx@smc.vnet.net
- Delivered-to: mathgroup-newsendx@smc.vnet.net
On Wed, Oct 02, 2013 @ 02:53 AM, Gebbo <nicolasgebbo at googlemail.com> 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@#}& @
Table[Last[worals@@(Delete[#,i]&/@{x,y,sx,sy})],{i,n}]
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 http://forums.wolfram.com/mathgroup/archive/2009/Feb/msg00478.html