滑动t检验

zjm1997

非原创，α=0.05。
需改动的部分见标黄部分，如研究时间段是1971~2016年，则N=46（时间跨度）、NYEAR=1971（起始年份）、K=10（步长，可根据自己需要设置K=5或K=10）。OPEN1是原始数据文件，OPEN2是计算结果会放到的文件。

α=0.05，步长为5时，其自由度为5+5-2=8，则t0.05=±2.31；步长为10时，其自由度为10+10-2=18，则t0.05=±2.10

PROGRAM N_MTT
PARAMETER (N=46 ,K=10 ,NYEAR=1971)
REAL *8 X(N),X1(N),X2(N),T(K:N-K+1)
REAL *8 V1,V2,EQ1,EQ2
OPEN (1,FILE="f:/spei01.txt")
DO I=1,N
READ (1,*) X(I)
ENDDO
CLOSE(1)

IF (K.EQ.5) THEN
A=2.31
ELSE IF (K.EQ.10) THEN
A=2.10
ELSE IF (K.EQ.12) THEN
A=2.05
ENDIF
DO L=K,N-K+1
DO I=1,K
X1(I)=X(L-K+I)
X2(I)=X(L+I-1)
ENDDO
CALL SVAR(X1,V1,K,N)
CALL SVAR(X2,V2,K,N)
CALL SAVAG(X1,EQ1,K,N)
CALL SAVAG(X2,EQ2,K,N)
T(L)=(EQ1-EQ2)/SQRT((2.0/K)*((K-1)*(V1+V2)/(K+K-2)))
PRINT *,NYEAR+L-1
ENDDO
OPEN (2,FILE="f:/mtt10spei01.txt")
DO L=K,N-K+1
WRITE (2,'(I4,3F18.2)') NYEAR+L-1,T(L),A,-A
ENDDO
CLOSE(2)
END
SUBROUTINE SVAR(X,V,K,N)
REAL *8 X(K),V
SUM=0
DO I=1,K
SUM=SUM+X(I)
ENDDO
AVAG=SUM/K
VAR=0
DO I=1,K
VAR=VAR+(X(I)-AVAG)**2
ENDDO
V=VAR/K
END
SUBROUTINE SAVAG(X,EQ,K,N)
REAL *8 X(K),EQ
SUMX=0
DO I=1,K
SUMX=SUMX+X(I)
ENDDO
EQ=SUMX/K
END

原始数据格式（就是只有一列）及计算结果

weidu699

为什么结果没有显示全部年份的数据

小白要逆袭

很棒，做论文正需要，谢谢楼主。

多洛洛

请教一下这个结果怎么分析呢

zjm1997

多洛洛 发表于 2020-4-24 21:30
请教一下这个结果怎么分析呢

可以在家园论坛里或者“鸠摩搜书”那个网站中搜索《现代气候统计诊断与预测技术(第二版）》这本书。

多洛洛

zjm1997 发表于 2020-4-25 07:22
可以在家园论坛里或者“鸠摩搜书”那个网站中搜索《现代气候统计诊断与预测技术(第二版）》这本书。

好的，买了这本书

omak718

非常感谢楼主~~~

pengtingting

T(L)=(EQ1-EQ2)/SQRT((2.0/K)*((K-1)*(V1+V2)/(K+K-2)))这里等号后面需要加一个负号嘛，目前跑出来的变大的年份突变点对应的是负值？