Hi,
Looking through the code for DERIVSIG, I'm a little confused about its
derivation.
From DERIV, it uses a 3-point Lagrange interpolant, so (as is in the
code), for inner points (ignoring end points for now)
;df/dx = y0*(2x-x1-x2)/(x01*x02)+y1*(2x-x0-x2)/(x10*x12)+y2*(2x-x0-x1 )/
(x20*x21)
; Where: x01 = x0-x1, x02 = x0-x2, x12 = x1-x2, etc.
In IDL, this is how it is written (from DERIV line 69-71,74-76):
x12 = x - shift(x,-1) ;x1 - x2
x01 = shift(x,1) - x ;x0 - x1
x02 = shift(x,1) - shift(x,-1) ;x0 - x2
d = shift(y,1) * (x12 / (x01*x02)) + $ ;Middle points
y * (1./x12 - 1./x01) - $
shift(y,-1) * (x01 / (x02 * x12))
If I do a standard error propagation formula
(sigma^2=sum(sigma_i^2*(\frac{\partial{f},\partial{i}),i) ), then I
get (in the terms defined above)
sigma_dfdx = sqrt( (2x-x1-x2)^2/(x01*x02)^2*sig_y0^2 + (2x-x0-x2)^2/
(x10*x12)^2*sig_y1^2 + (2x-x0-x1)^2/(x20*x21)^2*sig_y2^2)
In IDL, making use that x=x1, and since all the xij terms are squared
in the error propagation formula, we set xij=xji, I would get:
x01=shift(x,1)-x
x02=shift(x,1)-shift(x,-1)
x12=x-shift(x,-1)
sig_y0sqr=shift(sigy,1)*shift(sigy,1)
sig_y1sqr=sigy*sigy
sig_y2sqr=shift(sigy,-1)*shift(sigy,-1)
sigd=(x12^2/(x01*x02)^2*sig_y0sqr + (2*x-shift(x,-1)-shift(x,1))^2/
(x01*x12)^2*sig_y1sqr + x01^2/(x02*x12)^2*sig_y2sqr)
Yet in DERIVSIG, this is the equation used:
dsq=shift(x,-1)-shift(x,1)
dsq=dsq*dsq
dy=shift(y,-1)-shift(y,1)
sigd=(shift(sigy,-1)*shift(sigy,-1) + shift(sigy,1)*shift(sigy,1))/dsq
sigd=sqrt(sigd)
In the notation I used above, this is equivalent to saying:
sigd = sqrt( 1./x20^2*sig_y0^2 + 1./x20^2*sig_y2^2)
Working through how this compares with the equation I have above, the
only way I could get them to match is if DERVSIG is ASSUMING that x1-
x0 = x2-x1, or that the spacing is uniform (not unity, but equal on
each side of a point). This appears to be incorrect (the documentation
certainly doesn't mention that the spacing needs to be uniform). Have
I missed something somewhere?
Michael
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
