该子函数根据Dr. Jianping Li的fortran函数改编
function [fh,fl,dh,dl,dhl,tn,ta,ta95] = hcfx(n,x,f,coefh,coefl)
% For a time series f(n), computing:
% (1) mean in the high index years of x(n)
% (2) mean in the low index years of x(n)
% (3) composite difference between the mean of f(n) in high
% index years of x(n) and its climatology average
% (4) composite difference between the mean of f(n) in low
% index years of x(n) and its climate average
% (5) composite difference between the means of f(n) in high
% and low index years of x(n)
% input: n,x(n),f(n),coefh,coefl
% n: the length of time series
% x: control variable (index),eg：降水指数
% f: given series，eg:风场
% coefh: control parameter for high index (i.e., high index years are
% those in which x(i) > coefh)
% coefl: control parameter for low index (i.e., low index years are
% those in which x(i) < coefl)
% output: fh,fl,dh,dl,dhl,tn(5)
% fh: the mean of f in high index years of x(n)
% fl: the mean of f in low index years of x(n)
% dh: composite difference between the mean of f in high index years of x(n)
% and its climate mean (i.e., high index years minus cliamte mean).
% dl: composite difference between the mean of f in low index years of x(n)
% and its climate mean (i.e., low index years minus cliamte mean).
% dhl: composite difference between the means of f in high and low index years
% of x(n) (i.e., high minus low index years)
% tn(i,j): tn only equals 2., -2., 1., -1. or 0. corresponding to significant difference or not.
% tn=2. indicates that the difference is positive and significant.
% tn=-2. indicates that the difference is negative and significant.
% tn=1. indicates that the difference is positive but not significant.
% tn=-1. indicates that the difference is negative but not significant.
% tn=0. indicates the difference is zero.
% tn(1,j)~tn(5,j) are corresponding to the 90%,95%,98%,99% and 99.9% confident levels.
% j=1: siginificant test for dh
% j=2: siginificant test for dl
% j=3: siginificant test for dhl
% Feburary 11, 2002 by Jianping Li.
nh=length(f(x>=coefh));
hn=f(x>=coefh);
nl=length(f(x<=coefl));
ln=f(x<=coefl);
fh=mean(f(x>=coefh));
fl=mean(f(x<=coefl));
avef=mean(f);
dh=fh-avef;
dl=fl-avef;
dhl=fh-fl;
[tn(:,1),ta(:,1),ta95(:,1)]=diff_t_test(nh,n,hn,f);
[tn(:,2),ta(:,2),ta95(:,2)]=diff_t_test(nl,n,ln,f);
[tn(:,3),ta(:,3),ta95(:,3)]=diff_t_test(nh,nl,hn,ln);
end
function [tn,ta,ta95] = diff_t_test(n,m,x,y)
% n,m两个样本的样本数，n特殊年，m所有年
% x,y两个样本的序列
nn=10000;
tn(1:5)=0;
ax=mean(x);
sx=std(x);
vx=sx.^2;
ay=mean(y);
sy=std(y);
vy=sy.^2;
sxy=(nvx+mvy)/(n+m-2);
sxy=sxy*(1./n+1./m);
sxy=sqrt(sxy);
% if sxy==0
dxy=ax-ay;
if(dxy>0)
sn=2;
tn(1:5)=1;
else
sn=-2;
tn(1:5)=-1;
end
ta=abs(dxy)/sxy;
nm=n+m-2;
[ft(:,1),ft(:,2),ft(:,3),ft(:,4),ft(:,5)]=t_table(nn);
ta95=ft(nm,2);
for i=1:5
if(ta>=ft(nm,i))
tn(i)=sn;
end
end
end
function [ft90,ft95,ft98,ft99,ft999]=t_table(n)
% t-distribution, i.e., student's distribution
% t table with two-tailed (right and left tails) probabilities
% P(|t|>=ta)=a
% where a (alpha) significance level and (1-a)*100% is confindence level.
% t90: t-distribution test at 90% confidence level.
% t95: t-distribution test at 95% confidence level.
% t98: t-distribution test at 98% confidence level.
% t99: t-distribution test at 99% confidence level.
% t999: t-distribution test at 99.9% confidence level.
% By Dr. Jianping Li, January 5, 2000.
% n=10000;
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];
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];
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];
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];
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];
ft90=n90/1000;
ft95=n95/1000;
ft98=n98/1000;
ft99=n99/1000;
ft999=n999/1000;
for 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;
end
for 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;
end
for 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;
end
for 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;
end
end
网友评论