两段滑动t检验代码比较

HJY@E

网上找到了两段滑动t检验代码，但是同样的数据，结果有点点差异，如图1所示。

两种代码差距

图1

根据公式研究了一下，发现造成这种差异原因的是两段代码使用的方差计算方法不一样，如图2

两种方差方法

图2

哪种方法对呢，欢迎讨论^^

下附代码：

                               
登录/注册后可看大图

1. %% 滑动T检验，两种方法代码有删改
   
2. clc,clear;
   
3. %% 读取数据  
   
4. A=[1900 15.4;1901 14.6;1902 15.8;1903 14.8;1904 15;1905 15.1;1906 15.1;1907 15;1908 15.2;
   
5. 1909 15.4;1910 14.8;1911 15;1912 15.1;1913 14.7;1914 16;1915 15.7;1916 15.4;1917 14.5;
   
6. 1918 15.1;1919 15.3;1920 15.5;1921 15.1;1922 15.6;1923 15.1;1924 15.1;1925 14.9;1926 15.5;
   
7. 1927 15.3;1928 15.3;1929 15.4;1930 15.7;1931 15.2;1932 15.5;1933 15.5;1934 15.6;1935 15.1;
   
8. 1936 15.1;1937 16;1938 16;1939 16.8;1940 16.2;1941 16.2;1942 16;1943 15.6;1944 15.9;
   
9. 1945 16.2;1946 16.7;1947 15.8;1948 16.2;1949 15.9;1950 15.8;1951 15.5;1952 15.9;1953 16.8;
   
10. 1954 15.5;1955 15.8;1956 15;1957 14.9;1958 15.3;1959 16;1960 16.1;1961 16.5;1962 15.5;
    
11. 1963 15.6;1964 16.1;1965 15.6;1966 16;1967 15.4;1968 15.5;1969 15.2;1970 15.4;1971 15.6;
    
12. 1972 15.1;1973 15.8;1974 15.5;1975 16;1976 15.2;1977 15.8;1978 16.2;1979 16.2;1980 15.2;
    
13. 1981 15.7;1982 16;1983 16;1984 15.7;1985 15.9;1986 15.7;1987 16.7;1988 15.3;1989 16.1];
    
14. %% 方法1   来源https://blog.csdn.net/qq_44913577/article/details/106134923
    
15. myData=A;
    
16. timeSeries=myData(:,1);  % 时间序列数据
    
17. dataSeries=myData(:,2);  %  
    
18. dataCount=length(dataSeries);%检验数据长度
    
19. %% 设置步长与检验值
    
20. step=5;  % 步长  
    
21. v=step+step-2; % 计算自由度  
    
22. len1=step;%前一组样本数
    
23. len2=step; %后一组样本数
    
24. x=timeSeries(step:dataCount-step);
    
25. for i=step:dataCount-step   %设置每个基准点i   
    
26.     n1=dataSeries(i-step+1:i);   
    
27.     n2=dataSeries(i+1:i+step);     
    
28.     mean1=mean(n1);     
    
29.     mean2=mean(n2);      
    
30.     c=(len1+len2)/(len1*len2);      
    
31.     var1=1/len1*sum((n1-mean1).^2);     
    
32.     var2=1/len2*sum((n2-mean2).^2);     
    
33.     delta1=len1*var1+len2*var2;     
    
34.     delta=delta1/(len1+len2-2);      
    
35.     t1(i-step+1)=-(mean1-mean2)/sqrt(delta*c);
    
36. end
    
37. %% 制图
    
38. figure(1);
    
39. p1=plot(x,t1,'b-','linewidth',1.5,'DisplayName','方法1');  
    
40. xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);
    
41. ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
    
42. axis([min(x),max(x),min(t1),max(t1)]);
    
43. hold on  
    
44. 
    
45. %%方法2  来源http://bbs.06climate.com/forum.php?mod=viewthread&tid=49034&highlight=%BB%AC%B6%AFT%BC%EC%D1%E9
    
46. year=A(:,1);%year
    
47. y=A(:,2); %data
    
48. years=length(year);%数据的时长
    
49. %MMT方法检验%
    
50. %先设定子序列的长度
    
51. sublen=5;
    
52. n1=sublen;n2=sublen;
    
53. m=sqrt((1.0/n1)+(1.0/n2));
    
54. %开始计算每个基准点的统计量序列t
    
55. for i=sublen:years-sublen
    
56.     x1avg=mean(y(i-sublen+1:i));
    
57.     x2avg=mean(y(i+1:i+sublen));
    
58.     x1var=var(y(i-sublen+1:i));  
    
59.     x2var=var(y(i+1:i+sublen));
    
60.     s=sqrt((n1*x1var+n2*x2var)/(n1+n2-2));
    
61.     t(i)=(x2avg-x1avg)/(s*m);
    
62. end
    
63. i=sublen:years-sublen;
    
64. p2=plot(year(i),t(i),'r-','DisplayName','方法2');

复制代码

wjy_ecnu

都对。。。。。

HJY@E

wjy_ecnu 发表于 2022-5-4 15:17
都对。。。。。

个人倾向于小样本