R.G. Stockwell <noemail@please.com> wrote:
>> I am processing some long term time series of various data.
>> I computed the long term trends using the REGRESS fuction.
>> However, I need to also give the significance level of the
>> trends. I have been looking around for a code to do it but
>> have not found anything.
> Perhaps the sigma values will give you what you want. I don't remember
> offhand, but 1-sigma is the 67% significance level (or is that 90%?)
Here is what I use. Sorry for the messy pro... I was young then.
If you don't want to wade through the details, the important bit is:
; Compute the (2-tailed) t-statistic appropriate to the probability
t=t_cvf(prob/2.,reg.n/adj-2)
print,'A confidence interval for the slope is: ',reg.b+[-1,1]*reg.seb*t
(seb is standard-error-b, as I recall).
NOTE that if your data are auto-correlated then you will get the wrong
result unless you take this into account
-W.
; Explain the components of a regression generated by pp_regress
;
; Input
; reg - anonymous structure generated by pp_regress
; Options
; prob=x - probability value for 2-tailed test. Default is 0.05 (95% conf)
; /quiet - suppress info
; /simple - don't try to account for autocorr if it exists
; Output
; sig - return value, 1 if slope, b, is sig. diff from 0; 0 otherwise
; - Also writes info to the screen
; ac - value of first five auto correlations at lags 1...5
; acs - ac-sig: index of sig autocorrs, empty if there are none
; sig_ac - threshold used for significance
;
; A reasonable test of this routine is:
; n=24 & ntest=1000
; a=0 & x=indgen(n)
; for i=1,ntest do a=reg_explain(pp_regress(x,randomn(seed,n),/arr),prob=.1,/qu )+a
; print,a/float(ntest)
; Which should produce a number something like 0.10
;
; Nerdy example of use:
; print,'reg slope '+(['is not','is'])(reg_explain(/qu,pp_regress(x,y,/arr)))+' diff from zero'
function reg_explain1,reg,prob=prob,quiet=quiet,ac,acs,sig_ac
@comm_error
; Set up confidence level
if (not keyword_Set(prob)) then begin
prob=0.05
if (not keyword_Set(quiet)) then message,/cont,'Setting prob to 0.05 (2-tailed), ie 95% confidence'
endif
***some auto-corr stuff removed***
adj=1.0
; Compute the (2-tailed) t-statistic appropriate to the probability
t=t_cvf(prob/2.,reg.n/adj-2)
if (not keyword_Set(quiet)) then begin
print,'The series had '+shtstr(reg.n)+' elements & the x-range is: ',reg.xr
print,'The regression equation is: y='+shtstr(reg.b)+' * x + '+shtstr(reg.a)
print,'The r and r^2 value is: ',shtstr(reg.r)+', '+shtstr(reg.r^2)
print,'A confidence interval for the slope is: ',reg.b+[-1,1]*reg.seb*t
print,'The t-value used was: ',shtstr(t)
endif
; Decide if B is sig diff from zero, and return 1 if so else 0
if (abs(reg.b) gt reg.seb*t) then begin
if (not keyword_Set(quiet)) then print,'B is sig. diff. from zero'
return,1
endif else begin
if (not keyword_Set(quiet)) then print,'B is not sig. diff. from zero'
return,0
endelse
end
; ------------------------
function reg_explain,regs,ac,acs,sig_ac,_extra=e
for i=0,n_elements(regs)-1 do $
if (i eq 0) then ret=replicate(reg_explain1(regs(i),ac,acs,sig_ac,_extra=e),n _elements(regs)) $
else $
ret(i)=reg_explain1(regs(i),ac,acs,sig_ac,_extra=e)
if (n_elements(ret) eq 1) then ret=ret(0)
return,ret
end
--
William M Connolley | wmc@bas.ac.uk | http://www.antarctica.ac.uk/met/wmc/
Climate Modeller, British Antarctic Survey | Disclaimer: I speak for myself
I'm a .signature virus! copy me into your .signature file & help me spread!
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
