| Re: xerr [message #64353 is a reply to message #64276] |
Thu, 18 December 2008 11:51   |
laxsri
Messages: 6
Registered: November 2008
|
Junior Member |
|
|
Thanks all!
I shall look up the orthogonal regression Wox...I guess, I thought
this would have been implemented in IDL by somebody! (yea...lazy too)!
Good day
Lakshmi
On Dec 18, 9:35 pm, Wox <s...@nomail.com> wrote:
> On Wed, 17 Dec 2008 12:27:59 -0800 (PST), lakshmi <lax...@gmail.com>
> wrote:
>
>> I've been using mpfitfun to fit measured values of period (y) and
>> distances (x) in a linear equation y = a + bx.
>> I would like to know if we can include the measured uncertainties in x
>> values too?
>
> What you're looking for is called "Total least squares" or "orthogonal
> regression". Here is a reference + code you can use, translate to IDL
> and hopefully share it with use :-).
>
> Article:http://www.iop.org/EJ/abstract/0957-0233/18/11/025
> Matlab (off course) code:http://www.mathworks.com/matlabcentral/fileexchange/174 66
>
> It uses the method of the Lagrange Multipliers, so it should be
> possible to add constraints to the rico and intercept.
>
> If you're lazy: here's a simple code (not tested and without
> calculating errors) using info fromhttp://mathforum.org/library/drmath/view/63765.html. Change
> fixintercept to 1 when you want to fix the intercept.
>
> pro odr
> n=100
> x=findgen(n)+RANDOMN(seed,n)
> rico=1.2
> intercept=3
> fixintercept=0b
> y=rico*x+intercept+RANDOMN(seed,n)
>
> print,'Rico:',rico
> print,'Intercept:',intercept
>
> ; Centroid: orthogonal distance
> ; regression line goes through it
> n=n_elements(x)
> data=transpose([[x],[y]])
> centroid=total(data,2)/n
>
> ; Optional: Fix intercept
> if fixintercept then centroid=[0,intercept]
>
> data[0,*]-=centroid[0]
> data[1,*]-=centroid[1]
>
> SVDC, data, W, U, V
>
> smallest_singularvalue=min(W,ind)
> normal=reform(V[ind,*])
>
> rico=-normal[0]/normal[1]
> intercept=-rico*centroid[0]+centroid[1]
>
> print,'ODR...'
> print,'Rico:',rico
> print,'Intercept:',intercept
>
> window
> plot,x,y,psym=1
> oplot,x,rico*x+intercept
> end
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
