-
Notifications
You must be signed in to change notification settings - Fork 1
/
bootstraptrend.ncl
51 lines (43 loc) · 1.37 KB
/
bootstraptrend.ncl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
; This function is written to use bootstrap method to test the significance of trend
; Y[*] is the parameter you want to test, time[*] is the time series, n[1] is the times of repeating calculating.
; Written by Kai Zhang
function bootstraptrend(Y[*],time[*],n[1])
local time,sloptmp,b0,Ytmp,bstorage,pcdf,bounds,findind,zpdf,k,genindice,i,j,ntime
begin
ntime = dimsizes(time)
sloptmp = new((/ntime-1,ntime-1/),float)
Y1 = dtrend(Y,False)
;b0 = median(xi-xj/i-j)
do i = 0,ntime-2
do j = 1,ntime-1
sloptmp(i,j-1) = (Y1(j)-Y1(i))/(time(j)-time(i))
end do
end do
b0 = dim_median(ndtooned(sloptmp))
bstorage = new((/n/),float)
;generate bootstrap samples
do k = 0,n-1
sloptmp = 0
genindice = toint(random_uniform(0,ntime-1,ntime))
Ytmp = Y1(genindice)
do i = 0,ntime-2
do j = 1,ntime-1
sloptmp(i,j-1) = (Ytmp(j)-Ytmp(i))/(time(j)-time(i))
end do
end do
bstorage(k)= dim_median(ndtooned(sloptmp))
end do
opt = True
zpdf = pdfx(bstorage,100,opt) ;100 bins
pcdf = cumsum(zpdf,0)
bounds = zpdf@bin_bounds
findind = ind(bounds.ge.b0)
if (dimsizes(findind).gt.1) then
return(pcdf(findind(0)))
else
return(pcdf(0))
end if
end