滑动T检验

wjxbnd

不知道怎么出不了结果，这是我使用的程序：
C        THIS IS A PROGRAM FOR DETECTING ABRUPT CLIMATIC CHANGE
C        BY USING MOVING t-TEST TECHNIQUE
        PROGRAM MTT
        DIMENSION Y(1000),YN1(1000),YN2(1000),YC1(1000),YC2(1000),
     &        S(1000),T(1000),NY1(1000),NY2(1000),NNY1(1000),NNY2(1000)
        WRITE(*,10)
  10        FORMAT(5X,'N=70,IH=5,NYEAR=1930')
        READ(*,*)N,IH,NYEAR
C        **************************************************
C        * N:    SAMPLE SIZE                              *
C        * IH:   LENGTH OF SUB-SERIES                         *
C        * NYEAR: FIRST YEAR OF THE TIME SERIES           *
C        * Y(N):  ORIGINAL TIME SERIES                    *
C        * ************************************************
        OPEN(2,FILE='d:\2.DAT')
        READ(2,*)(Y(I),I=1,N)
        N1=N-IH+1
        N2=N-2*IH+1
        C1=0.0
        C2=0.0
        DO 20 I=1,IH
  20        C1=C1+Y(I)
        DO 30 I=IH+1,2*IH
  30        C2=C2+Y(I)
        DO 40 I=1,N-IH
        D1=C1-Y(I)+Y(I+IH)
        YN1(I)=C1/IH
  40        C1=D1
        DO 50 I=IH+1,N-IH
        D2=C2-Y(I)+Y(I+IH)
        YN2(I-IH)=C2/IH
  50        C2=D2
        YN1(N1)=C1/IH
        YN2(N2)=C2/IH
        DO 60 I=1,N1
        YC1(I)=0.0
        DO 70 J=I,IH+I-1
  70        YC1(I)=YC1(I)+(Y(J)-YN1(I))*(Y(J)-YN1(I))
  60        CONTINUE
        DO 80 I=1,N2
        YC2(I)=0.0
        DO 90 J=IH+I,2*IH+I-1
  90        YC2(I)=YC2(I)+(Y(J)-YN2(I))*(Y(J)-YN2(I))
  80        CONTINUE
        DO 100 I=1,N2
  100        S(I)=SQRT((YC1(I)+YC2(I))/(IH+IH-2))
        DO 110 I=1,N2
  110        T(I)=(YN1(I)-YN2(I))/(S(I)*SQRT(2.0/IH))
        DO 120 I=1,N2
        NY1(I)=NYEAR+I-1
        NNY1(I)=NY1(I)+IH-1
        NY2(I)=NYEAR+I+IH-1
  120        NNY2(I)=NY2(I)+IH-1
        IF(IH.EQ.5)THEN
        A=4.0
        B=-4.0
        ELSE IF(IH.EQ.10)THEN
        A=3.2
        B=-3.2
        ELSE
        A=3.0
        B=-3.0
        END IF
        OPEN(3,FILE='d:\3.DAT',STATUS='NEW')
        WRITE(3,130)
  130        FORMAT(30X,'ABRUPT CLIMATIC CHANGE ANALYSIS'/)
        WRITE(3,140)
  140        FORMAT(25X,'TM',4X,'t-TEST(0.01)'/)
        DO 150 I=1,N2
        WRITE(3,160)NY1(I),NNY1(I),NY2(I),NNY2(I),T(I),A,B
  160        FORMAT(1X,I4,'-',I4,'--',I4,'-',I4,1X,3F8.2)
  150        CONTINUE
        OPEN(4,FILE='d:\4.DAT',STATUS='NEW')
        DO 170 I=1,N2
        WRITE(4,180)NNY1(I),T(I),A,B
  180        FORMAT(1X,I4,3F8.2)
  170        CONTINUE
        STOP
        END
问题如图
还请高手指点！

river

错误提示说的很明白了，就是说你要生成的那个文件已经存在了，但是OPEN(3,FILE='d:\3.DAT',STATUS='NEW')红色的部分要求它不能被覆盖，所以您老人家把,STATUS='NEW'去掉好了

wjxbnd

你真是我的救星，可以用了。非常感谢！

wjxbnd

大侠，给留个QQ吧。以后要多向你学习！

Legendary

楼主没有觉得这个程序的自由度错了？当子序列长度为5的时候，自由度应该是5+5-2=8，查表得+-3.36，而程序里是自由度为+-4

chuanhao

wjxbnd 发表于 2014-4-5 09:43
你真是我的救星，可以用了。非常感谢！

你这个程序，我运行不出来呢？？

wjxbnd

程序一般不会有问题，只要数据路径及设置正确是可以的啊。

coco4046

Legendary 发表于 2014-9-8 23:08
楼主没有觉得这个程序的自由度错了？当子序列长度为5的时候，自由度应该是5+5-2=8，查表得+-3.36，而程序里 ...

请问自由度是怎么算的？我一直没搞明白，原始序列-2 ？

墨荷

自由度=n1+n2-2=5+5-2=8

想请教下楼主，为什么有明显增大的转折点是负值

libiecho

学习了