| Re: xerr [message #64365 is a reply to message #64276] |
Thu, 18 December 2008 02:35  |
Wout De Nolf
Messages: 194
Registered: October 2008
|
Senior Member |
|
|
On Wed, 17 Dec 2008 12:27:59 -0800 (PST), lakshmi <laxsri@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/17466
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 from
http://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
