- 积分
- 1935
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2012-5-8
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
求f(n)在指数x(n)为高指数年(x(n)>coefh的年）的平均值fh、低指数年(x(n)<coefh的年）的平均值fl、高指数年与气候平均的合成差dh、低指数年与气候平均的合成差dl、以及高低指数年的合成差dhl和差的显著性tn(5,3)。
c-----*----------------------------------------------------6---------7--
c For a time series f(n), computing:
c (1) mean in the high index years of x(n)
c (2) mean in the low index years of x(n)
c (3) composite difference between the mean of f(n) in high
c index years of x(n) and its climatology average
c (4) composite difference between the mean of f(n) in low
c index years of x(n) and its climate average
c (5) composite difference between the means of f(n) in high
c and low index years of x(n)
c input: n,x(n),f(n),coefh,coefl
c n: the length of time series
c x: control variable (index)
c f: given series
c coefh: control parameter for high index (i.e., high index years are
c those in which x(i) > coefh)
c coefl: control parameter for low index (i.e., low index years are
c those in which x(i) < coefl)
c output: fh,fl,dh,dl,dhl,tn(5)
c fh: the mean of f in high index years of x(n)
c fl: the mean of f in low index years of x(n)
c dh: composite difference between the mean of f in high index years of x(n)
c and its climate mean (i.e., high index years minus cliamte mean).
c dl: composite difference between the mean of f in low index years of x(n)
c and its climate mean (i.e., low index years minus cliamte mean).
c dhl: composite difference between the means of f in high and low index years
c of x(n) (i.e., high minus low index years)
c tn(i,j): tn only equals 2., -2., 1., -1. or 0. corresponding to significant difference or not.
c tn=2. indicates that the difference is positive and significant.
c tn=-2. indicates that the difference is negative and significant.
c tn=1. indicates that the difference is positive but not significant.
c tn=-1. indicates that the difference is negative but not significant.
c tn=0. indicates the difference is zero.
c tn(1,j)~tn(5,j) are corresponding to the 90%,95%,98%,99% and 99.9% confident levels.
c j=1: siginificant test for dh
c j=2: siginificant test for dl
c j=3: siginificant test for dhl
c Feburary 11, 2002 by Jianping Li.
subroutine differencehl1(n,x,f,coefh,coefl,fh,fl,dh,dl,dhl,tn)
dimension x(n),f(n),tn(5,3)
real hn(n),ln(n) !work array
call meanvar(n,f,avef,sf,vf)
fh=0.
fl=0.
nh=0
nl=0
do k=1,n
if(x(k).gt.coefh)then
nh=nh+1
hn(nh)=f(k)
endif
if(x(k).lt.coefl)then
nl=nl+1
ln(nl)=f(k)
endif
enddo
call meanvar1(nh,hn,fh)
call meanvar1(nl,ln,fl)
dh=fh-avef
dl=fl-avef
dhl=fh-fl
call diff_t_test(nh,n,hn,f,tn(1,1))
call diff_t_test(nl,n,ln,f,tn(1,2))
call diff_t_test(nh,nl,hn,ln,tn(1,3))
return
end
c-----*----------------------------------------------------6---------7--
c Student's t-test for the difference between two means
c input: n,m,x(n),y(m)
c n: number of x
c m: number of y
c x(n): series
c y(m): series
c output: tn(5)
c tn: tn only equals 2, -2., 1, -1. or 0. corresponding to significant difference or not.
c tn=2. indicates that the difference is positive and significant.
c tn=-2. indicates that the difference is negative and significant.
c tn=1. indicates that the difference is positive but not significant.
c tn=-1. indicates that the difference is negative but not significant.
c tn=0. indicates the difference is zero.
c tn(1)~tn(5) are corresponding to the 90%,95%,98%,99% and 99.9% confident levels.
c By Dr. Jianping Li, March 05, 2004
subroutine diff_t_test(n,m,x,y,tn)
parameter(nn=10000)
real x(n),y(m),tn(5)
real ft(nn,5) !work array
do i=1,5
tn(i)=0.
enddo
call meanvar(n,x,ax,sx,vx)
call meanvar(m,y,ay,sy,vy)
sxy=(n*vx+m*vy)/(n+m-2.)
sxy=sxy*(1./n+1./m)
sxy=sqrt(sxy)
if(sxy.eq.0.)return
dxy=ax-ay
if(dxy.gt.0.)then
sn=2.
do i=1,5
tn(i)=1.
enddo
else
sn=-2.
do i=1,5
tn(i)=-1.
enddo
endif
ta=abs(dxy)/sxy
nm=n+m-2
call t_table(ft(1,1),ft(1,2),ft(1,3),ft(1,4),ft(1,5))
do i=1,5
if(ta.ge.ft(nm,i))then
tn(i)=sn
endif
enddo
return
end
c-----*----------------------------------------------------6---------7--
c Computing the mean ax, standard deviation sx
c and variance vx of a series x(i) (i=1,...,n).
c input: n and x(n)
c n: number of raw series
c x(n): raw series
c output: ax, sx and vx
c ax: the mean value of x(n)
c sx: the standard deviation of x(n)
c vx: the variance of x(n)
c By Dr. LI Jianping, May 6, 1998.
subroutine meanvar(n,x,ax,sx,vx)
dimension x(n)
ax=0.
vx=0.
sx=0.
do 10 i=1,n
ax=ax+x(i)
10 continue
ax=ax/float(n)
do 20 i=1,n
vx=vx+(x(i)-ax)**2
20 continue
vx=vx/float(n)
sx=sqrt(vx)
return
end
c-----*----------------------------------------------------6---------7--
c Computing the mean ax of a series x(i) (i=1,...,n).
c input: n and x(n)
c n: number of raw series
c x(n): raw series
c output: ax
c ax: the mean value of x(n)
c By Dr. LI Jianping, May 6, 1999.
subroutine meanvar1(n,x,ax)
dimension x(n)
ax=0.
vx=0.
sx=0.
do 10 i=1,n
ax=ax+x(i)
10 continue
ax=ax/float(n)
return
end
c-----*----------------------------------------------------6---------7--
c t-distribution, i.e., student's distribution
c t table with two-tailed (right and left tails) probabilities
c P(|t|>=ta)=a
c where a (alpha) significance level and (1-a)*100% is confindence level.
c t90: t-distribution test at 90% confidence level.
c t95: t-distribution test at 95% confidence level.
c t98: t-distribution test at 98% confidence level.
c t99: t-distribution test at 99% confidence level.
c t999: t-distribution test at 99.9% confidence level.
c By Dr. Jianping Li, January 5, 2000.
subroutine t_table(ft90,ft95,ft98,ft99,ft999)
parameter(n=10000)
dimension ft90(n),ft95(n),ft98(n),ft99(n),ft999(n)
dimension n90(30),n95(30),n98(30),n99(30),n999(30)
data n90/ 6314,2920,2353,2132,2015,1943,1895,1860,1833,1812,
* 1796,1782,1771,1761,1753,1746,1740,1734,1729,1725,
* 1721,1717,1714,1711,1708,1706,1703,1701,1699,1697/
data n95/12706,4303,3182,2776,2571,2447,2365,2306,2262,2228,
* 2201,2179,2160,2145,2131,2120,2110,2101,2093,2086,
* 2080,2074,2069,2064,2060,2056,2052,2048,2045,2042/
data n98/31821,6965,4541,3747,3365,3143,2998,2896,2821,2764,
* 2718,2681,2650,2624,2602,2583,2567,2552,2539,2528,
* 2518,2508,2500,2492,2485,2479,2473,2467,2462,2457/
data n99/63657,9925,5841,4604,4032,3707,3499,3355,3250,3169,
* 3106,3055,3012,2977,2947,2921,2898,2878,2861,2845,
* 2831,2819,2807,2797,2787,2779,2771,2763,2756,2750/
data n999/636619,31598,12941,8610,6859,5959,5405,5041,4781,
* 4587,4437,4318,4221,4140,4073,4015,3965,3922,3883,3850,
* 3819,3792,3767,3745,3725,3707,3690,3674,3659,3646/
do 10 i=1,30
ft90(i)=n90(i)/1000.
ft95(i)=n95(i)/1000.
ft98(i)=n98(i)/1000.
ft99(i)=n99(i)/1000.
ft999(i)=n999(i)/1000.
10 continue
do 20 i=31,40
fi=(i-30.)/10.
ft90(i)=1.697-(1.697-1.684)*fi
ft95(i)=2.042-(2.042-2.021)*fi
ft98(i)=2.457-(2.457-2.423)*fi
ft99(i)=2.750-(2.750-2.704)*fi
ft999(i)=3.646-(3.646-3.551)*fi
20 continue
do 30 i=41,60
fi=(i-40.)/20.
ft90(i)=1.684-(1.684-1.671)*fi
ft95(i)=2.021-(2.021-2.000)*fi
ft98(i)=2.423-(2.423-2.390)*fi
ft99(i)=2.704-(2.704-2.660)*fi
ft999(i)=3.551-(3.551-3.460)*fi
30 continue
do 40 i=61,120
fi=(i-60.)/60.
ft90(i)=1.671-(1.671-1.658)*fi
ft95(i)=2.000-(2.000-1.980)*fi
ft98(i)=2.390-(2.390-2.358)*fi
ft99(i)=2.660-(2.660-2.617)*fi
ft999(i)=3.460-(3.460-3.373)*fi
40 continue
do 50 i=121,n
fi=(i-120.)/(n-120.)
ft90(i)=1.658-(1.658-1.645)*fi
ft95(i)=1.980-(1.980-1.960)*fi
ft98(i)=2.358-(2.358-2.326)*fi
ft99(i)=2.617-(2.617-2.576)*fi
ft999(i)=3.373-(3.373-2.291)*fi
50 continue
return
end
|
|
窥视卡
雷达卡
发表于 2012-12-1 20:53:33
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
