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