| Theil–Sen estimator anyone? [message #82896] |
Sat, 19 January 2013 23:41  |
bstecklu
Messages: 14
Registered: February 2012
|
Junior Member |
|
|
Is there an IDL implementation of the Theil–Sen estimator, also known as Kendall
robust line-fit method? I'm aware of robust_poly_fit.pro but that's not quite
the same.
Bringfried
|
|
|
|
| Re: Theil-Sen estimator anyone? [message #87470 is a reply to message #82896] |
Wed, 05 February 2014 05:37  |
bstecklu
Messages: 14
Registered: February 2012
|
Junior Member |
|
|
On Sunday, January 20, 2013 8:41:26 AM UTC+1, bstecklu wrote:
> Is there an IDL implementation of the Theil-Sen estimator, also known as Kendall
>
> robust line-fit method? I'm aware of robust_poly_fit.pro but that's not quite
>
> the same.
>
>
>
> Bringfried
To answer my own question here is a basic implementation.
; perform Theil-Sen robust linear regression
; cf. https://en.wikipedia.org/wiki/Theil%E2%80%93Sen_estimator#CI TEREFSen1968
function ts,x,y
n=n_elements(x)
xx=x#replicate(1,n)
xx=xx-transpose(xx)
yy=y#replicate(1,n)
yy=yy-transpose(yy)
ts=yy/(xx+(xx eq 0))
idx=where(xx ne 0)
ts=ts[idx]
;sort
ts=ts[sort(ts)]
; slope
a=median(ts)
; constant
b=median(y)-a*median(x)
np=n*(n-1)/2.
; 95% confidence interval for slope, cf. http://pubs.usgs.gov/tm/2006/tm4a7/
z=1.96*sqrt(n*(n-1)*(2*n+5)/18.)
ru=ts[round((np+z)/2.+1)]
rl=ts[round((np-z)/2.)]
;done
return,[a,b,rl,ru]
end
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
