C------------------------------------------------------------------------ C C Program name: DIST C C Description: For integer-valued censored Normal(mu,sigma) data, C given estimates of mu and sigma, determine the C exact distribution of X-bar and R for particular C sample sizes. C C Last revised: July 1989 (currently under development) C C----------------------------------------------------------------------- C REAL*8 MU, SIGMA, PHI(11), + VAL1(10), P1(10), SUM1, + VAL2(19), P2(19), SUM2, + VAL3(28), P3(28), SUM3, + VAL4(37), P4(37), SUM4, + VAL5(46), P5(46), SUM5, + VAL6(55), P6(55), SUM6, + P(10), ONE, F(10), S(10,10), + M(10,10), PROBR(10), SUMP INTEGER NUMI, MINI, N C WRITE(*,10) 10 FORMAT(////,15X,'Distributions of Xbar and R',//, + 15X,' from rounded data',/////) C 22 WRITE(*,20) 20 FORMAT(//,1X,'Enter the estimates of mu and ', + 'sigma and return ',/ + ,1X,'----> ') READ(*,*) MU, SIGMA C IF (SIGMA.LE.1.1D0) GOTO 23 WRITE(*,*) ' Sigma must be no more than 1.1' CALL PAUSE GOTO 22 23 CALL GETPHI(NUMI,MINI,MU,SIGMA,PHI) CALL GETCDF(NUMI,PHI,P,ONE,F) CALL GETSIJ(NUMI,F,S) C CALL GETP1(NUMI,PHI,VAL1,P1,SUM1) CALL GETP2(NUMI,P1,VAL2,P2,SUM2) CALL GETP3(NUMI,P1,P2,VAL3,P3,SUM3) CALL GETP4(NUMI,P2,VAL4,P4,SUM4) CALL GETP5(NUMI,P2,P3,VAL5,P5,SUM5) CALL GETP6(NUMI,P3,VAL6,P6,SUM6) C 25 WRITE(*,30) 30 FORMAT(//,1X,'Enter the sample size (up to 6) and return ',/ + ,1X,'(enter 0 to quit)',/ + ,1X,'----> ') READ(*,*) N C IF (N.LE.6) GOTO 100 WRITE(*,*) ' Sample size must be no more than 6' CALL PAUSE GOTO 25 C 100 GOTO (110,120,130,140,150,160),N GOTO 999 C 110 CALL OUTX(NUMI,MINI,MU,SIGMA,VAL1,P1,SUM1,N) GOTO 200 120 CALL OUTX(NUMI,MINI,MU,SIGMA,VAL2,P2,SUM2,N) GOTO 200 130 CALL OUTX(NUMI,MINI,MU,SIGMA,VAL3,P3,SUM3,N) GOTO 200 140 CALL OUTX(NUMI,MINI,MU,SIGMA,VAL4,P4,SUM4,N) GOTO 200 150 CALL OUTX(NUMI,MINI,MU,SIGMA,VAL5,P5,SUM5,N) GOTO 200 160 CALL OUTX(NUMI,MINI,MU,SIGMA,VAL6,P6,SUM6,N) C 200 CALL GETMIJ(NUMI,S,N,M) CALL GETDSN(NUMI,P,N,M,PROBR,SUMP) CALL OUTR(NUMI,N,PROBR,SUMP) C GOTO 25 C 999 STOP END C C---------------------------------------------------------------------- C SUBROUTINE GETPHI(NUMI,MINI,MU,SIGMA,PHI) C REAL*8 MU, SIGMA, SHIFTMU, NORMLO, NORMHI, PHI(11), Z, NORMCP INTEGER MINI, MAXI, NUMI, NUMIP1 C NORMLO = -4.0D0*SIGMA + MU + 0.5D0 NORMHI = 4.0D0*SIGMA + MU - 0.5D0 C IF (NORMLO.LT.0.0D0.AND.NORMHI.LT.0.0D0) THEN MINI = DINT(NORMLO-1.0D0) MAXI = DINT(NORMHI) ELSE IF (NORMLO.LT.0.0D0.AND.NORMHI.GE.0.0D0) THEN MINI = DINT(NORMLO-1.0D0) MAXI = DINT(NORMHI+1.0D0) ELSE IF (NORMLO.GE.0.0D0.AND.NORMHI.GE.0.0D0) THEN MINI = DINT(NORMLO) MAXI = DINT(NORMHI+1.0D0) ELSE WRITE(*,*) ' ERROR IN FINDING NUMI' END IF C NUMI = MAXI - MINI + 1 NUMIP1 = NUMI + 1 SHIFTMU = MU + (1.0D0 - DBLE(MINI)) C DO 100 I=1,NUMIP1 Z = ((I - 0.5D0) - SHIFTMU) / SIGMA PHI(I) = NORMCP(Z) 100 CONTINUE C RETURN END C C-------------------------------------------------------------------- C DOUBLE PRECISION FUNCTION NORMCP(Z) C C This function returns the standard normal cumulative C probability of being less than Z. C REAL*8 Z,C(7),Y,Q,S DATA C/.319381530D0,-.356563782D0,1.781477937D0, + -1.821255978D0,1.330274429D0,.2316419D0, + 2.506628725D0/ IF (Z .GT. 6.0D0) THEN NORMCP = 1.0D0 C ELSEIF (Z .LT. -6.0D0) THEN NORMCP = 0.0D0 C ELSE Y=Z IF (Z .LT. 0.0D0) Y=-Z Q=1.0D0/(1.0D0+C(6)*Y) S=((((C(5)*Q+C(4))*Q+C(3))*Q+C(2))*Q+C(1))*Q NORMCP = S*DEXP(-Y*Y/2.0D0)/C(7) IF (Z .GT. 0.0D0) NORMCP=1.0D0-NORMCP ENDIF C 100 RETURN END C C----------------------------------------------------- C C SUBROUTINE GETCDF(NUMI,PHI,P,ONE,F) C REAL*8 PHI(11), P(10), ONE, F(10) INTEGER NUMI C P(1) = PHI(2) - PHI(1) F(1) = P(1) C ONE = P(1) DO 10 I=2,NUMI P(I) = PHI(I+1) - PHI(I) F(I) = F(I-1) + P(I) ONE = ONE + P(I) 10 CONTINUE C RETURN END C C C---------------------------------------------------- C SUBROUTINE GETSIJ(NUMI,F,S) C REAL*8 F(10), S(10,10) INTEGER NUMI C DO 20 I=1,NUMI DO 10 J=1,NUMI S(I,J) = 0.0D0 10 CONTINUE 20 CONTINUE C DO 25 J=1,NUMI S(1,J) = F(J) 25 CONTINUE C DO 40 I=2,NUMI DO 30 J=I,NUMI S(I,J) = F(J) - F(I-1) 30 CONTINUE 40 CONTINUE C RETURN END C C----------------------------------------------------- C SUBROUTINE GETP1(NUMI,PHI,VAL1,P1,SUM1) C REAL*8 PHI(25), VAL1(10), P1(10), SUM1 INTEGER NUMI C VAL1(1) = 1 P1(1) = PHI(2) - PHI(1) C SUM1 = P1(1) DO 10 I=2,NUMI VAL1(I) = DBLE(I) P1(I) = PHI(I+1) - PHI(I) SUM1 = SUM1 + P1(I) 10 CONTINUE C RETURN END C C C---------------------------------------------------- C SUBROUTINE GETP2(NUMI,P1,VAL2,P2,SUM2) C REAL*8 P1(10), VAL2(19), P2(19), SUM2 INTEGER NUMI, NUMIP1, MAX2, LAST C SUM2 = 0.0D0 NUMIP1 = NUMI + 1 MAX2 = (2*NUMI) - 1 C DO 50 L=1,NUMI VAL2(L) = (L+1.0D0)/2.0D0 P2(L) = 0.0D0 DO 40 I=1,L P2(L) = P2(L) + ( P1(I) * P1(L-I+1) ) 40 CONTINUE SUM2 = SUM2 + P2(L) 50 CONTINUE C DO 70 L=NUMIP1,MAX2 VAL2(L) = (L+1.0D0)/2.0D0 P2(L) = 0.0D0 LAST = MAX2 - L + 1 DO 60 I=1,LAST P2(L) = P2(L) + ( P1(NUMI-LAST+I) * P1(NUMI-I+1) ) 60 CONTINUE SUM2 = SUM2 + P2(L) 70 CONTINUE C RETURN END C C C---------------------------------------------------- C SUBROUTINE GETP3(NUMI,P1,P2,VAL3,P3,SUM3) C REAL*8 P1(10), P2(19), VAL3(28), P3(28), SUM3 INTEGER NUMI, MAX2, MAX2P1, MAX3, LAST C SUM3 = 0.0D0 MAX2 = (2*NUMI) - 1 MAX2P1 = MAX2 + 1 MAX3 = (3*NUMI) - 2 C DO 50 L=1,MAX2 VAL3(L) = (L+2.0D0)/3.0D0 P3(L) = 0.0D0 LAST = MIN0(NUMI,L) DO 40 I=1,LAST P3(L) = P3(L) + ( P1(I) * P2( L-I+1 ) ) 40 CONTINUE SUM3 = SUM3 + P3(L) 50 CONTINUE C DO 70 L=MAX2P1,MAX3 VAL3(L) = (L+2.0D0)/3.0D0 P3(L) = 0.0D0 LAST = MAX3 - L + 1 DO 60 I=1,LAST P3(L) = P3(L) + ( P1( NUMI-LAST+I ) * P2( MAX2-I+1 ) ) 60 CONTINUE SUM3 = SUM3 + P3(L) 70 CONTINUE C RETURN END C C---------------------------------------------------- C SUBROUTINE GETP4(NUMI,P2,VAL4,P4,SUM4) C REAL*8 P2(19), VAL4(37), P4(37), SUM4 INTEGER NUMI, MAX2, MAX2P1, MAX4 C SUM4 = 0.0D0 MAX2 = (2*NUMI) - 1 MAX2P1 = MAX2 + 1 MAX4 = (4*NUMI) - 3 C DO 50 L=1,MAX2 VAL4(L) = (L+3.0D0)/4.0D0 P4(L) = 0.0D0 DO 40 I=1,L P4(L) = P4(L) + ( P2(I) * P2(L-I+1) ) 40 CONTINUE SUM4 = SUM4 + P4(L) 50 CONTINUE C DO 70 L=MAX2P1,MAX4 VAL4(L) = (L+3.0D0)/4.0D0 P4(L) = 0.0D0 LAST = MAX4 - L + 1 DO 60 I=1,LAST P4(L) = P4(L) + ( P2(MAX2-LAST+I) * P2(MAX2-I+1) ) 60 CONTINUE SUM4 = SUM4 + P4(L) 70 CONTINUE C RETURN END C C---------------------------------------------------- C SUBROUTINE GETP5(NUMI,P2,P3,VAL5,P5,SUM5) C REAL*8 P2(19), P3(28), VAL5(46), P5(46), SUM5 INTEGER NUMI, MAX2, MAX3, MAX3P1, MAX5, LAST C SUM5 = 0.0D0 MAX2 = (2*NUMI) - 1 MAX3 = (3*NUMI) - 2 MAX3P1 = MAX3 + 1 MAX5 = (5*NUMI) - 4 C DO 50 L=1,MAX3 VAL5(L) = (L+4.0D0)/5.0D0 P5(L) = 0.0D0 LAST = MIN0(MAX2,L) DO 40 I=1,LAST P5(L) = P5(L) + ( P2(I) * P3( L-I+1 ) ) 40 CONTINUE SUM5 = SUM5 + P5(L) 50 CONTINUE C DO 70 L=MAX3P1,MAX5 VAL5(L) = (L+4.0D0)/5.0D0 P5(L) = 0.0D0 LAST = MAX5 - L + 1 DO 60 I=1,LAST P5(L) = P5(L) + ( P2( MAX2-LAST+I ) * P3( MAX3-I+1 ) ) 60 CONTINUE SUM5 = SUM5 + P5(L) 70 CONTINUE C RETURN END C C---------------------------------------------------- C SUBROUTINE GETP6(NUMI,P3,VAL6,P6,SUM6) C REAL*8 P3(28), VAL6(55), P6(55), SUM6 INTEGER NUMI, MAX3, MAX3P1, MAX6 C SUM6 = 0.0D0 MAX3 = (3*NUMI) - 2 MAX3P1 = MAX3 + 1 MAX6 = (6*NUMI) - 5 C DO 50 L=1,MAX3 VAL6(L) = (L+5.0D0)/6.0D0 P6(L) = 0.0D0 DO 40 I=1,L P6(L) = P6(L) + ( P3(I) * P3(L-I+1) ) 40 CONTINUE SUM6 = SUM6 + P6(L) 50 CONTINUE C DO 70 L=MAX3P1,MAX6 VAL6(L) = (L+5.0D0)/6.0D0 P6(L) = 0.0D0 LAST = MAX6 - L + 1 DO 60 I=1,LAST P6(L) = P6(L) + ( P3(MAX3-LAST+I) * P3(MAX3-I+1) ) 60 CONTINUE SUM6 = SUM6 + P6(L) 70 CONTINUE C RETURN END C C---------------------------------------------------- C SUBROUTINE GETMIJ(NUMI,S,N,M) C REAL*8 S(10,10), SN(10,10), M(10,10) INTEGER NUMI, N C DO 20 I=1,NUMI DO 10 J=1,NUMI SN(I,J) = S(I,J)**N M(I,J) = 0.0D0 10 CONTINUE 20 CONTINUE C DO 40 I=1,NUMI DO 30 J=I+1,NUMI M(I,J) = SN(I,J) - SN(I+1,J) + - SN(I,J-1) + SN(I+1,J-1) 30 CONTINUE 40 CONTINUE C RETURN END C C---------------------------------------------------- C SUBROUTINE GETDSN(NUMI,P,N,M,PROBR,SUMP) C REAL*8 P(10), M(10,10), PROBR(10), SUMP INTEGER NUMI, N, NUMIMR C PROBR(1) = 0.0D0 DO 10 I=1,NUMI PROBR(1) = PROBR(1) + P(I)**N 10 CONTINUE C SUMP = PROBR(1) DO 40 L=2,NUMI PROBR(L) = 0.0D0 NUMIMR = NUMI - (L-1) DO 20 I=1,NUMIMR PROBR(L) = PROBR(L) + M(I,I+(L-1)) 20 CONTINUE SUMP = SUMP + PROBR(L) 40 CONTINUE C RETURN END C C---------------------------------------------------- C SUBROUTINE OUTX(NUMI,MINI,MU,SIGMA,VALX,PX,SUMX,N) C INTEGER NUMI, MINI, N, MAXX REAL*8 MU, SIGMA, VAL, VALX(55), PX(55), SUMX C WRITE(*,4) 4 FORMAT(5(/)) WRITE(*,5) ' MU: ', MU WRITE(*,5) ' SIGMA: ', SIGMA 5 FORMAT(1X,A8,F12.6) CALL PAUSE C WRITE(*,30) N 30 FORMAT(5(/),1X,'Distribution of Xbar (N=',I2,'): ',// +10X,'Xbar Value',10X,'Probability',//) MAXX = (N*NUMI) - (N-1) DO 40 I=1,MAXX VAL = MINI + VALX(I) - 1.0D0 WRITE(*,20) VAL, PX(I) 40 CONTINUE WRITE(*,28) SUMX CALL PAUSE C 20 FORMAT(10X,F10.5,10X,F15.12) 28 FORMAT(30X,'---------------',/ +30X,F15.12) C RETURN END C C------------------------------------------------------- C SUBROUTINE OUTR(NUMI,N,PROBR,SUMP) C REAL*8 PROBR(10), SUMP INTEGER NUMI, N, R C WRITE(*,80) N 80 FORMAT(5(/),1X,'Distribution of the Range ( N = ',I2,'): ',// +10X,'R Value ',10X,'Probability',//) DO 85 I=1,NUMI R=I-1 WRITE(*,90) R, PROBR(I) 90 FORMAT(16X,I4,10X,F15.12) 85 CONTINUE WRITE(*,68) SUMP 68 FORMAT(30X,'---------------',/ +30X,F15.12) CALL PAUSE C RETURN END C C------------------------------------------------------- C SUBROUTINE PAUSE C CHARACTER*1 GO C C Pause to let the user view a display. C WRITE(*,70) 70 FORMAT(/,1X,'STRIKE CARRIAGE RETURN ', +'WHEN YOU WISH TO CONTINUE ') READ(*,75) GO 75 FORMAT(A1) C RETURN END C C---------------------------------------------------- C SUBROUTINE CHK(NUM) C INTEGER NUM C WRITE(*,*) 'CHECKPOINT ',NUM C RETURN END C C----------------------------------------------------