home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / sigm / vol269 / regrs1.for < prev    next >
Encoding:
Text File  |  1986-05-22  |  32.1 KB  |  682 lines
  1. C [REGRS1.FOR]
  2. C
  3. C     ** MULTIPLE LINEAR REGRESSION ANALYSIS WITH INTERVAL ESTIMATION **
  4. C        (╛▌╣▓ ╝▐¡│╢▓╖ ╠▐▌╛╖ ─ ╕╢▌╜▓├▓)
  5. C
  6. C    Yoshio MONMA (JUG CP/M No.43)            85-03-21    
  7. C    1514-63 Takata-cho, Khohoku-ku
  8. C    Yokohama 223, Japan
  9. C
  10. C
  11.     PROGRAM    REGRS1
  12. C
  13. C       The following declaration are for the machines with 56k TPA.
  14. C    You can increase the data and variable sizes as
  15. C       '80' ---> '100' and '13' ---> '16'
  16. C    if your system has 64k TPA.
  17.  
  18.       INTEGER*1    BEL,SI
  19.       INTEGER*2    ITRS(13),JIN(13),RP80
  20.       REAL*4       CMNT(18),FMT(18),VARS(18)
  21.       REAL*4       X(80,13),S(13,13),R(13,13),XX(20,13)    
  22.       REAL*4       AVGX(13),SIGX(13),TOL(13),YY(20)             
  23.       REAL*4       YEST(80),RES(80),WRK(80),XMIN(13),XMAX(13)
  24.       REAL*8       DAY,SAMP
  25. C                                                                       
  26.       EQUIVALENCE  (ITRS(1),JIN(1)), (RES(1),WRK(1))
  27. C
  28.       DATA         BEL,RP80,SI/Z'07',Z'1B40',Z'0F'/, LX,LS/80,13/
  29.       data day/'85-03-21'/
  30. C
  31.       WRITE(1,1000)
  32.  1000   FORMAT(1H1,5X,'** Multiple Linear Regression Analysis **'/)
  33.    10 WRITE(1,1010)
  34.  1010   FORMAT(1H ,8X,'Specify input file (FORT0X.DAT, X in I1) = ')
  35.       READ(1,1020) INP0
  36.  1020   FORMAT(I1)
  37.       IF (INP0.LT.6.OR.INP0.GT.9)     GOTO 10
  38. C    Set RP-80 to condensed mode
  39.       WRITE(2,2000) SI
  40.  2000   FORMAT(1H ,A1)
  41. C            * Get initial parameters (14I5)
  42.       READ(INP0,500) IN,JPY,MPL,INP1,KDTA,ISWP,IRES,JGRH,INTV,
  43.      1               ICFL1,IPRB1,ICFL2,IPRB2,NPRD                     
  44.   500   FORMAT(16I5)
  45. C
  46. C    You may modify '80' --> '100' and '12' ---> '15' in the following
  47. C       statement, if your system has 64k TPA.
  48.       IF (IN.GT.80.OR.JPY.GT.12.OR.MPL.GT.5) GOTO 199
  49.       IF (INTV.GT.5.OR.NPRD.GT.20) GOTO 199   
  50. C                If the system has the feature
  51. C     CALL DATE(DAY)                                                    
  52. C                                                                       
  53.       READ(INP0,510) CMNT,SAMP
  54.   510   FORMAT(1X,18A4,A8)
  55. C                                                                       
  56.       WRITE(2,2010) DAY,SAMP
  57.  2010   FORMAT(1H /11X,'** Multiple Linear Regression Analysis with ',
  58.      1  'Interval Estimation **'//14X,'Yoshio MONMA',60X,'Date ',A8/14X,
  59.      2  'JUG CP/M No.43'//11X,
  60.      3  'Program Name = REGRS1',30X,'Sample = ',A8)
  61.       CALL BIGCHR(85,SAMP,2)
  62.       WRITE(1,1030) SAMP
  63.  1030   FORMAT(1H ,50X,'Sample = ',A8)
  64.       WRITE(1,1040) CMNT
  65.  1040   FORMAT(1H ,5X,18A4)
  66.       WRITE(2,1040) CMNT
  67.       WRITE(1,1050) IN,JPY,INP1,KDTA,ISWP,IRES,JGRH,
  68.      1              INTV,ICFL1,IPRB1,ICFL2,IPRB2,NPRD                    
  69.  1050   FORMAT(1H0,'* Initial Parameters *'/11X,'IN =',I4,
  70.      1    ', JPY =',I3,', INP1 =',I2,', KDTA =',I2,', ISWP =',I2,    
  71.      2    ', IRES =',I2,', JGRH =',I2/11X,'INTV =',I2,', ICFL1 =',I3,   
  72.      3    ', IPRB1 =',I3,', ICFL2 =',I3,', IPRB2 =',I3,', NPRD =',I3/)
  73.       WRITE(2,1050) IN,JPY,INP1,KDTA,ISWP,IRES,JGRH,
  74.      1              INTV,ICFL1,IPRB1,ICFL2,IPRB2,NPRD                    
  75. C
  76.       WRITE(2,1060) 
  77.  1060   FORMAT(1H0,10X,'* Description of Variables *'/7X,
  78.      1           'Var No',10X, 'Description')                                  
  79.       DO 20 J=1,JPY                                                     
  80.          READ(INP0,510) VARS
  81.          WRITE(2,1070) J,VARS
  82.  1070      FORMAT(1H ,I10,10X,18A4)
  83.    20 CONTINUE                                                          
  84.       JP = JPY-1
  85.       JPY1 = JPY+1                                                      
  86.       IF (MPL.GT.1) JPY1 = JPY1+MPL-1                                   
  87.       IF (INP0.NE.INP1) REWIND INP1
  88.       READ (INP0,510) FMT
  89.       DO 30 I=1,IN                                                      
  90.          READ(INP1,FMT) (X(I,J),J=1,JPY)                  
  91.          X(I,JPY1) = 1.0                                                
  92.    30 CONTINUE                                                          
  93.       IF (NPRD.EQ.0)                    GOTO 35                         
  94. C            * Get the data for prediction
  95.       READ(INP0,510) FMT
  96.       READ(INP0,FMT) ((XX(I,J),J=1,JP),I=1,NPRD)
  97. C            * Transformation of variables
  98.    35 READ(INP0,500) (ITRS(J),J=1,JPY)
  99.       DO 44 J=1,JPY                                                     
  100.          DO 40 I=1,IN                                                   
  101.             WRK(I) = X(I,J)                                            
  102.    40    CONTINUE                                                       
  103.          CALL TRNSV0(INP0,IN,WRK,ITRS(J))
  104.          DO 42 I=1,IN                                                   
  105.             X(I,J) = WRK(I)                                            
  106.    42    CONTINUE                                                       
  107.    44 CONTINUE                                                          
  108.       IF (NPRD.EQ.0)                    GOTO 50                         
  109.       DO 50 J=1,JP                                                      
  110.          DO 45 I=1,IN                                                   
  111.             WRK(I) = XX(I,J)                                           
  112.    45    CONTINUE                                                       
  113.          CALL TRNSV0(INP0,NPRD,WRK,ITRS(J))
  114.          DO 50 I=1,IN                                                   
  115.             XX(I,J) = WRK(I)                                           
  116.    50 CONTINUE                                                          
  117. C                                                                       
  118.       IF (MPL.LE.1)                     GOTO 60                         
  119.       WRITE(2,1080) MPL
  120.  1080   FORMAT(1H0,10X,'* Polynomial Regression, Order =',I2,' *')
  121.       JPY = JPY+MPL-1                                                   
  122.       MPL = JPY-1                                                       
  123.       DO 55 I=1,IN                                                      
  124.          X(I,JPY) = X(I,2)                                              
  125.          DO 53 J=2,MPL                                                  
  126.             X(I,J) = X(I,1)**J                                          
  127.    53    CONTINUE                                                       
  128.    55 CONTINUE                                                          
  129.       IF (NPRD.EQ.0)                    GOTO 60                         
  130.       DO 59 I=1,IN                                                      
  131.          DO 57 J=2,MPL                                                  
  132.             XX(I,J) = XX(I,1)**J                                        
  133.    57    CONTINUE                                                       
  134.    59 CONTINUE                                                          
  135. C                                                                       
  136.    60 IF (KDTA.NE.1) CALL MATPRI('Data Matrix     ',X,1,LX,IN,JPY)              
  137. C                                                                       
  138. C                * Copy the data matrix to FORT10.DAT
  139.       REWIND 10
  140.       DO 66 I=1,IN                                                      
  141.          DO 63 J=1,JPY                                                  
  142.             WRK(J) = X(I,J)                                             
  143.    63    CONTINUE                                                       
  144.          WRITE(10) (WRK(J),J=1,JPY)                    
  145.    66 CONTINUE                                                          
  146. C
  147. C                    * SS matrix (═▓╬│▄Ñ╛╖▄ ╖▐«│┌┬)
  148.       CALL SSSP(IN,JPY,S,LS,WRK)
  149. C
  150.       DO 70 J=1,JPY1                                                    
  151.          TOL(J) = 0.00001*S(J,J)                                        
  152.    70 CONTINUE                                                          
  153. C                * Statistics
  154.       DF = IN-1
  155.       DO 75 J=1,JPY 
  156.          AVGX(J) = S(JPY1,J)                                            
  157.          SIGX(J) = SQRT(S(J,J)/DF)                                      
  158.          CALL MINMAX(X(1,J),IN,XMIN(J),XMAX(J))                         
  159.    75 CONTINUE                                                          
  160. C
  161.       CALL VECPRI('Mean    ',AVGX,1,JPY)
  162.       CALL VECPRI('Sigma   ',SIGX,1,JPY)
  163.       CALL VECPRI('Min.    ',XMIN,1,JPY)
  164.       CALL VECPRI('Max.    ',XMAX,1,JPY)
  165.       IF (NPRD.GT.0) CALL MATPRI('For Prediction  ',XX,1,20,NPRD,JP)
  166.       IF (ISWP.EQ.1) CALL MATPRI('S.S. Matrix     ',S,1,LS,JPY1,JPY1)   
  167. C                                                                       
  168. C                    * Correlation matrix
  169.       CALL CORREL(S,LS,JPY,R,LS)                                        
  170.       CALL MATPRI('Correl. Matrix  ',R,1,LS,JPY,JPY)             
  171. C                                                                       
  172.       IF (ISWP.EQ.1) WRITE(2,1090)
  173.  1090   FORMAT(1H0,10X,'* Sweep Process *')
  174.       WRITE(1,1090)
  175. C                                                                       
  176.       JP = JPY-1                                                        
  177.       SYY = S(JPY,JPY)                                                  
  178.       DO 80 J=1,JP                                                      
  179.          JIN(I) = -1                                                    
  180.    80 CONTINUE                                                          
  181.       JIN(JPY) = 0                                                      
  182.       JIN(JPY1) = 1                                                     
  183.       DFE = DF                                                          
  184.       DO 85 J=1,JP                                                      
  185.          IF (S(J,J).LT.TOL(J))          GOTO 85                         
  186.          CALL SWEEP1(S,LS,J,JPY1)                                       
  187.          JIN(J) = 1                                                     
  188.          DFE = DFE-1.0                                                  
  189.       IF (ISWP.EQ.1) CALL MATPRI('1/S             ',S,1,LS,JPY1,JPY1)
  190.    85 CONTINUE                                                          
  191. C                * Save data matrix and reg. coef.
  192.       WRITE(1,1095)
  193.  1095   FORMAT(1H0,5X,'* Saving ... data matrix & reg. coef.')
  194.       REWIND 10
  195.       WRITE(10) IN,JPY,IRES,JGRH
  196.       WRITE(10) CMNT,SAMP
  197.       DO 90 J=1,JPY
  198.         WRITE(10) (X(I,J),I=1,IN)
  199.    90 CONTINUE
  200.       WRITE(10) S(JPY1,JPY),(S(J,JPY),J=1,JP)
  201. C
  202. C                 * Analysis of variance
  203.       CALL ANOVA(S,LS,JIN,SIGX,IN,JPY,SYY,DFE,SE)                       
  204. C
  205.       IF (IRES.EQ.1) GOTO 100
  206.       WRITE(2,1100)
  207.  1100   FORMAT(1H1/)
  208.       CALL BIGCHR(85,SAMP,2)
  209.       WRITE(2,1040) CMNT
  210. C                    * Estimate
  211.       CALL ESTIM(S,LS,X,LX,JIN,IN,JPY,DFE,YEST,RES,SE,            
  212.      1           INTV,ICFL1,IPRB1,ICFL2,IPRB2,XX,NPRD,YY)              
  213. C                                                                       
  214. C
  215.       CALL DWRTIO(RES,IN,DWR)                                           
  216.       WRITE(2,2080) DWR
  217.  2080   FORMAT(1H0,10X,'* Durbin-Watson ratio *'//14X,'DWR =',F7.4)       
  218.   100 WRITE(2,1100)
  219.       STOP
  220. C
  221.   199 WRITE(1,1900) BEL
  222.  1900   FORMAT(1H ,A1,'*** Error in the Initial Parameter(s)!')
  223.       STOP ERROR
  224. C                                                                       
  225.       END
  226. C
  227.       SUBROUTINE   BIGCHR(VPOS,STRING,MODE)
  228. C
  229. C    ** Print STRING in Enlarged Mode at VPOS **
  230. C
  231. C    Written by Yoshio MONMA on 85-03-20
  232. C
  233. C     Arguments:
  234. C        VPOS      Vertical position
  235. C        STRING    Character string (A8), must be given in exact length
  236. C        MODE      Mode to be set after the printing STRING
  237. C             = 0  Standard (Pica, 10char/inch)
  238. C             = 1  Elite (Elite, 12char/inch) 
  239. C             = 2  Condensed (15char/inch)
  240. C
  241. C    This routine is for EPSON RP-80 Printer.
  242. C
  243.       INTEGER*1    BIGM,ESC,LETF,NUL,SI,SO
  244.       INTEGER*1    VPOS
  245.       REAL*8       STRING
  246. C
  247.       DATA         ESC/Z'1B'/,  NUL/'00'/, SI/Z'0F'/, SO/Z'0E'/
  248.       DATA         BIGM/Z'4D'/, LETF/Z'66'/
  249. C
  250.       WRITE(2,200) ESC,LETF,NUL,VPOS,SO,STRING
  251.   200   FORMAT(1H ,5A1,A8)
  252.       IF (MODE.EQ.1) WRITE(2,210) ESC,BIGM
  253.   210   FORMAT(1H ,2A1)
  254.       IF (MODE.EQ.2) WRITE(2,210) SI
  255.       RETURN
  256.       END
  257. C
  258.       SUBROUTINE   SSSP(II,JJ,S,LS,X)
  259. C
  260. C     ** Sum of Squares and Sum of Products **
  261. C
  262. C     * Reference, T.Haga & S.Hashimoto: ó╢▓╖╠▐▌╛╖ ─ ╝¡╛▓╠▐▌╠▐▌╛╖ú, P.27
  263. C
  264. C    Arguments:
  265. C    II    Data ize
  266. C    JJ    No. of variables
  267. C    S    S.S. matrix
  268. C    X    Work area
  269. C
  270. C    This routine reads the data matrix from FORT10.DAT.
  271. C
  272.       REAL*4       S(LS,1),X(1)
  273. C
  274.       JJ1 = JJ+1
  275.       DO 14 J=1,JJ
  276.          S(JJ1,J) = 0.0
  277.          DO 12 J1=1,J
  278.             S(J,J1) = 0.0
  279.    12    CONTINUE
  280.    14 CONTINUE
  281. C
  282.       REWIND 10
  283.       DO 20 I=1,II
  284.          READ(10) (X(J),J=1,JJ)
  285.          FI = I
  286.          F = (FI-1.0)/FI
  287.          DO 18 J=1,JJ
  288.             X(J) = X(J)-S(JJ1,J)
  289.             S(JJ1,J) = S(JJ1,J)+X(J)/FI
  290.             DO 16 J1=1,J
  291.                S(J,J1) = S(J,J1)+X(J)*X(J1)*F
  292.    16       CONTINUE
  293.    18    CONTINUE
  294.    20 CONTINUE
  295. C
  296.       DO 24 J=1,JJ
  297.          S(J,JJ1) = -S(JJ1,J)
  298.          DO 22 J1=1,J
  299.             S(J1,J) = S(J,J1)
  300.    22    CONTINUE
  301.    24 CONTINUE
  302. C
  303.       S(JJ1,JJ1) = 1.0/FI
  304. C
  305.       RETURN
  306.       END
  307. C
  308.       SUBROUTINE   ANOVA(S,LS,JIN,SIGX,IN,JPY,SYY,DFE,SE)               
  309. C                                                                       
  310. C     ** Analysis of Variance and Multiple Correlation **               
  311. C                                                                       
  312. C     * Reference, T.Haga & S.Hashimoto: ó╢▓╖╠▐▌╛╖ ─ ╝¡╛▓╠▐▌╠▐▌╛╖ú, P.27
  313. C
  314. C     * Arguments
  315. C     S         Sweeped out SS matrix
  316. C     LS        Size of S
  317. C     JIN       JIN(I)=1 if X(J) is adopted for a predictor variable,
  318. C               otherwise JIN(J)=-1                                  
  319. C     SIGX      Standard deviation
  320. C     IN        No. of sample (Data size)
  321. C     JPY       No. of variables (JP+1)                              
  322. C     SYY       Totak sum of squares for Y(I)                        
  323. C     DFE       Degree of freedom for residuals
  324. C     SE        Residual sum of squares
  325. C                                                                       
  326. C     * REFERENCE, T.Haga & S.Hashimoto: ╢▓╖╠▐▌╛╖ ─ ╝¡╛▓╠▐▌╠▐▌╛╖, P.48  
  327. C                                                                       
  328.       INTEGER*2    JIN(1)
  329.       REAL*4       S(LS,LS),SIGX(1)                             
  330. C                                                                       
  331.       FIN = IN                                                          
  332.       JP = JPY-1                                                        
  333.       JPY1 = JPY+1                                                      
  334.       SE = S(JPY,JPY)                                                   
  335.       SR  = SYY-SE                                                      
  336.       DF  = IN-1                                                        
  337.       DFR = DF-DFE                                                      
  338.       IF (DFR.EQ.0.0) DFR = 0.0001                                      
  339.       VR  = SR/DFR                                                      
  340.       VE  = SE/DFE                                                      
  341.       VT = SYY/DF                                                       
  342.       SEE = SQRT(VE)                                                    
  343.       F0 = VR/VE                                                        
  344. C                                                                       
  345.       IJ1 = IN-JP-1
  346.       F95 = PF(0.025,JP,IJ1)
  347.       F99 = PF(0.005,JP,IJ1)
  348.       JTST = '  '                                                       
  349.       IF (F0.GE.F95) JTST = '* '                                        
  350.       IF (F0.GE.F99) JTST = '**'                                        
  351.       WRITE(1,100) SR, DFR, VR, F0, JTST,                               
  352.      1             SE,DFE,VE,SEE,                                       
  353.      2             SYY,DF,VT                                            
  354.   100 FORMAT(1H0/11X,'* Analysis of Variance *'//6X,70('-') /10X,'S.V.',
  355.      1       7X,'S.S.',7X,'D.F.',5X,'M.S.',8X,'F0',5X,'F-test',3X,      
  356.      2       'Sig(E)'  /6X,70('-') /6X,'Regression',F11.4,F9.0,2F11.4,  
  357.      3       4X,A2/6X,'Residual',F13.4,F9.0,F11.4,15X,F12.4 /6X,70('-')
  358.      4      /6X,'Total',F16.4,F9.0,F11.4 /6X,70('-')/)                  
  359. C                                                                       
  360.          WRITE(2,100) SR, DFR, VR, F0, JTST,
  361.      1             SE, DFE, VE, SEE,                                       
  362.      2             SYY, DF, VT
  363.       R2 = 1.0-SE/SYY                                                   
  364.       RR2 = 1.0-VE/VT                                                   
  365.       RRR2 = 1.0-(2.0*FIN-DFE)/(2.0*FIN-DF)*(VE/VT)                     
  366.       R = SQRT(R2)                                                      
  367.       RR = SQRT(RR2)                                                    
  368.       RRR = SQRT(RRR2)                                                  
  369.       WRITE(1,110) R,R2,RR,RR2,RRR,RRR2                                 
  370.   110 FORMAT(1H0,10X,'* Multiple Correlation *'//31X,'R',9X,'R**2' /6X, 
  371.      1       'Ordinary',13X,F6.4,6X,F6.4 /6X,'Adjusted',13X,F6.4,6X,    
  372.      2       F6.4 /6X,'Double Adjusted',6X,F6.4,6X,F6.4/)               
  373. C                                                                       
  374.       WRITE(2,110) R,R2,RR,RR2,RRR,RRR2                                 
  375.       WRITE(1,120)                                                      
  376.   120 FORMAT(1H0,10X,'* Regression Coefficients and t-Test *'/13X,      
  377.      1 'Y = B(0) + B(1)*X(1) + B(2)*X(2) + ... + B(P)*X(P)'//10X,'J',   
  378.      2 7X,'B(J)',10X,'Std-B(J)',5X,'Sig(B)',5X,'t(B)',5X,'t-Test'/)     
  379.       WRITE(2,120)                                                      
  380.       T95 = PT(0.025,IJ1)
  381.       T99 = PT(0.005,IJ1)
  382.       DO 10 J=1,JP                                                      
  383.          IF (JIN(J).LT.0)               GOTO 10                         
  384.          B = S(J,JPY)                                                   
  385.          BB = B*SIGX(JPY)/SIGX(J)                                       
  386.          SIGB = SQRT(VE*S(J,J))                                         
  387.          T = B/SIGB                                                     
  388.          JTST = '  '                                                    
  389.          IF (ABS(T).GE.T95) JTST = '* '                                 
  390.          IF (ABS(T).GE.T99) JTST = '**'                                 
  391.          WRITE(1,130) J,B,BB,SIGB,T,JTST                                
  392.   130      FORMAT(1H ,I10,F14.5,F15.5,F11.4,F10.4,5X,A2)
  393.          WRITE(2,130) J,B,BB,SIGB,T,JTST                                
  394.    10 CONTINUE                                                          
  395.       IF (JIN(JPY1).GT.0) WRITE(1,140) S(JPY1,JPY)                      
  396.   140   FORMAT(1H0,'Constant =',F14.5)
  397.       IF (JIN(JPY1).GT.0) WRITE(2,140) S(JPY1,JPY)                      
  398.       RETURN                                                            
  399.       E N D                                                             
  400. C
  401.       SUBROUTINE  ESTIM(S,LS,X,LX,JIN,IN,JPY,DFE,YEST,RES,SE,     
  402.      1                  INTV,ICFL1,IPRB1,ICFL2,IPRB2,XX,NPRD,YY)       
  403. C                                                                       
  404. C     ** Estimate, Residuals and Interval Estimation **
  405. C
  406. C    Written by Yoshio MONMA (JUG-CP/M No.43)
  407. C
  408. C     * REFERENCE, T.Haga & S.Hashimoto: ó╢▓╖╠▐▌╛╖ ─ ╝¡╛▓╠▐▌╠▐▌╛╖ú, P.189
  409. C
  410. C     * Arguments                                                       
  411. C     S            Sweeped out SS matrix
  412. C     LS           Size of S
  413. C     X(LX,JPY)    Data matrix
  414. C     LX           Maximum data size
  415. C     JIN(J)       Index for variables picked up
  416. C     IN           Data size
  417. C     JPY          No. of variables (p+1)                               
  418. C     DFE          Degree of freedom for residuals
  419. C     RES(I)       Residuals                                             
  420. C     SE           Residual sum of squares
  421. C     INTV = 0     Average estimate only
  422. C          = 1     Avg. + CI (╝▌╫▓-╕╢▌)
  423. C          = 2     Avg. + SCI (─▐│╝▐-╝▌╫▓-╕╢▌)
  424. C          = 3     Avg. + PI (╓┐╕-╕╢▌)
  425. C          = 4     Avg. + TI (╖«╓│-╕╢▌)
  426. C          = 5     Avg. + STI (─▐│╝▐-╖«╓│-╕╢▌)
  427. C     ICFL        Confidence evel (%)
  428. C     IPRB        Probability level (%)
  429. C     XX(I,J)     Specified valies of predictor variables
  430. C     NPRD        No. of specified points (NPRD <= 20)
  431. C     YY(I)       Estimate for XX(I,J)
  432. C
  433.       INTEGER*2    JIN(1)
  434.       REAL*4       S(LS,1),X(LX,1),YEST(1),RES(1)                
  435.       REAL*4       XX(LS,1),YY(NPRD)                                    
  436.       REAL*8       RESTAB(5)
  437. C
  438.       DATA    RESTAB/'* Table ','of Resid','uals and',' Interva',
  439.      1                 'ls *    '/
  440. C
  441.       JPY1 = JPY+1                                                      
  442.       JP = JPY-1                                                        
  443.       NDF = DFE                                                         
  444.       T95 = PT(0.025,NDF)                                               
  445.       T99 = PT(0.005,NDF)                                               
  446.       VE = SE/DFE                                                       
  447.       ALPH1 = (100.0-ICFL1)/100.0                                       
  448.       ALPH2 = (100.0-ICFL2)/100.0                                       
  449.       GMA1  = (100.0-IPRB1)/100.0                                       
  450.       GMA2  = (100.0-IPRB2)/100.0                                       
  451. C                                                                       
  452.       INTV1 = INTV+1
  453.       GOTO (10,11,12,13,14,15), INTV1
  454. C
  455.    10 WRITE(1,110) RESTAB(1),RESTAB(2),RESTAB(5)
  456.   110   FORMAT(1H0,10X,2A8,'ua',A8)
  457.       WRITE(2,110) RESTAB(1),RESTAB(2),RESTAB(5)
  458.                     GOTO 20
  459.    11 WRITE(1,111) RESTAB
  460.   111 FORMAT(1H0,10X,3A8,' Confidence',2A8) 
  461.       WRITE(2,111) RESTAB
  462.                     GOTO 20
  463.    12 WRITE(1,112) RESTAB
  464.   112 FORMAT(1H0,10X,3A8,' Simultaneous Confidence',2A8)
  465.       WRITE(2,112) RESTAB
  466.                     GOTO 20
  467.    13 WRITE(1,113) RESTAB
  468.   113 FORMAT(1H0,10X,3A8,' Prediction',2A8)
  469.       WRITE(2,113) RESTAB
  470.                      GOTO 20
  471.    14 WRITE(1,114) RESTAB
  472.   114 FORMAT(1H0,10X,3A8,' Tolerance',2A8/)
  473.       WRITE(2,114) RESTAB
  474.                     GOTO 20
  475.    15 WRITE(1,115) RESTAB
  476.   115 FORMAT(1H0,10X,3A8,' Simultaneous Tolerance',2A8)
  477.       WRITE(2,115) RESTAB
  478. C
  479.    20 IF (INTV.GE.1.AND.INTV.LE.3) WRITE(2,220) ICFL1,ICFL2
  480.   220  FORMAT(1H0,13X,'Confidence Level =',I3,'% and',I3,'%'/)
  481.       IF (INTV.GE.4) WRITE(2,224) ICFL1,IPRB1,ICFL2,IPRB2
  482.   224   FORMAT(1H0,13X,'Confidence Level =',I3,'% with Probability of',
  483.      1       I3,'%')                                                    
  484.       IF (INTV.EQ.0) WRITE(2,230)
  485.   230   FORMAT(1H0,9X,'I',4X,'Y(obs)',4X,'Y(est)',3X,'Residual',2X,
  486.      1       'Sig(res)',5X,'t',6X,'t-Test'/)
  487.       IF (INTV.NE.0) WRITE(2,235) ICFL2,ICFL1,ICFL1,ICFL2               
  488.   235   FORMAT(1H0,9X,'I',4X,'Y(obs)',4X,'Y(est)',3X,'Residual',2X,
  489.      1       'Sig(res)',5X,'t',6X,'t-Test',12X,'Lower',14X,'Upper'/,    
  490.      2       76X,I3,'%',I10,'%',I8,'%',I10,'%'/)                        
  491. C                                                                       
  492.       DO 30 I=1,IN                                                      
  493.          YEST(I) = 0.0                                                  
  494. C     * C = 1/N + DM2/(N-1), DM2 = Mahalanobis' Generalized Distance
  495.          C = 0.0                                                        
  496.          DO 25 J=1,JPY1                                                 
  497.             IF (JIN(J).LE.0)            GOTO 25                         
  498.             YEST(I) = YEST(I) + X(I,J)*S(J,JPY)                         
  499.             DO 22 J1=1,JPY1                                             
  500.                IF (JIN(J1).LE.0)        GOTO 22                         
  501.                C = C + X(I,J)*S(J,J1)*X(I,J1)                           
  502.    22       CONTINUE                                                    
  503.    25    CONTINUE                                                       
  504.          RES(I) = X(I,JPY)-YEST(I)                                      
  505.          C1 = 1.0-C                                                     
  506.          SIGRES = SQRT(C1*(SE-RES(I)*RES(I)/C1)/(DFE-1.0))              
  507.          T = RES(I)/SIGRES                                              
  508.          JTST = '  '                                                    
  509.          IF (ABS(T).GT.T95) JTST = '* '                                 
  510.          IF (ABS(T).GT.T99) JTST = '**'                                 
  511. C                                       * Interval Estimation
  512.          CALL INTRVL(INTV,ALPH1,ALPH2,GMA1,GMA2,IN,NDF,C,VE,YEST(I),    
  513.      1               YL1,YL2,YU1,YU2)                              
  514.          IF (INTV.NE.0) WRITE(2,240) I,X(I,JPY),YEST(I),RES(I),SIGRES,  
  515.      1                               T,JTST,YL2,YL1,YU1,YU2             
  516.   240      FORMAT(1H ,I10,5F10.4,5X,A2,5X,4F10.4)                      
  517.          IF (INTV.EQ.0) WRITE(2,245) I,X(I,JPY),YEST(I),RES(I),SIGRES,  
  518.      1                               T,JTST                             
  519.   245    FORMAT(1H ,I10,5F10.4,5X,A2)
  520.          WRITE(10) YEST(I),RES(I),YL2,YL1,YU1,YU2   
  521.    30 CONTINUE                                                          
  522. C                                                                       
  523.       IF (NPRD.EQ.0)                    GOTO 50
  524.       WRITE(1,120)                                                      
  525.   120   FORMAT(1H0,10X,'* Prediction for specified X *'/)
  526.       WRITE(2,120)
  527.       DO 50 I=1,NPRD                                                    
  528.          YY(I) = 0.0                                                    
  529.          DO 35 J=1,JP                                                   
  530.             IF (JIN(J).LE.0)            GOTO 35                         
  531.             YY(I) = YY(I) + XX(I,J)*S(J,JPY)                            
  532.    35    CONTINUE                                                       
  533.          YY(I) = YY(I) + S(JPY1,JPY)                                    
  534.          XX(I,JPY) = YY(I)                                              
  535.          XX(I,JPY1) = 1.0                                               
  536.          C = 0.0                                                        
  537.          DO 45 J=1,JPY1                                                 
  538.                IF (JIN(J).LE.0)         GOTO 45                         
  539.             DO 40 J1=1,JPY1                                             
  540.                IF (JIN(J1).LE.0)        GOTO 40                         
  541.                C = C + XX(I,J)*S(J,J1)*XX(I,J1)                         
  542.    40       CONTINUE                                                    
  543.    45    CONTINUE                                                       
  544.          CALL INTRVL(INTV,ALPH1,ALPH2,GMA1,GMA2,IN,NDF,C,VE,YY(I),      
  545.      1               YL1,YL2,YU1,YU2)                              
  546.          IF (INTV.EQ.0) WRITE(2,250) I,YY(I)                            
  547.   250      FORMAT(1H ,I10,10X,F10.4)
  548.          IF (INTV.NE.0) WRITE(2,255) I,YY(I),YL2,YL1,YU1,YU2
  549.   255      FORMAT(1H ,I10,10X,F10.4,30X,12X,4F10.4)
  550.    50 CONTINUE                                                          
  551. C
  552.       RETURN                                                            
  553.       E N D                                                             
  554. C
  555.       SUBROUTINE   INTRVL(INTV,ALPH1,ALPH2,GMA1,GMA2,N,NDF,C,VE,YES,    
  556.      1                    YL1,YL2,YU1,YU2)                         
  557. C                                                                       
  558. C     ** Interval Estmation **                                         
  559. C                                                                       
  560. C    Written by Yoshio MONMA (JUG-CP/M No.43)
  561. C
  562. C     * Arguments                                                       
  563. C     INTV = 0: No Interval Estimation
  564. C          = 1: Confidence Interval
  565. C          = 2: Simultaneous Confidence Interval
  566. C          = 3: Prediction Interval
  567. C          = 4: Tolerance Interval
  568. C          = 5: Simultaneous Tolerance Interval
  569. C     ALPH1,ALPH2  Probability of t-Distribution (conf. level)          
  570. C     GMA1,GMA2    Probability of normal distribution
  571. C     N         Data size (no. of data points)
  572. C     NDF       Degree of freedom (N-P-1) 
  573. C     C       C = 1/n + DM2/(n-1), DM2 = Mahalanobis' Generalized Distance
  574. C     VE        Variance of residuals
  575. C     YES       Estimate
  576. C     YL1,YL2   Lower limits of the interval
  577. C     YU1,YU2   Upper limits of the interval
  578. C                                                                       
  579. C
  580.       IF (INTV.EQ.0) RETURN
  581.       GOTO (10,20,30,40,50), INTV                                       
  582. C                                                                       
  583.    10 AL12 = ALPH1/2.0
  584.       AL22 = ALPH2/2.0
  585.       PT12 = PT(AL12,NDF)
  586.       PT22 = PT(AL22,NDF)
  587.       X0  = SQRT(C*VE)                                                  
  588.    15 X1 = PT12*X0                                                      
  589.       X2 = PT22*X0                                                      
  590.                                         GOTO 60                         
  591.    20 NNDF = N-NDF
  592.       PF1 = PF(ALPH1,NNDF,NDF)
  593.       PF2 = PF(ALPH2,NNDF,NDF)
  594.       X0  = SQRT(C*VE)                                                  
  595.       X1  = SQRT((N-NDF)*PF1)*X0                                        
  596.       X2  = SQRT((N-NDF)*PF2)*X0                                        
  597.                                         GOTO 60                         
  598.    30 AL12 = ALPH1/2.0
  599.       AL22 = ALPH2/2.0
  600.       PT12 = PT(AL12,NDF)
  601.       PT22 = PT(AL22,NDF)
  602.       X0 = SQRT((C+1.0)*VE)                                             
  603.                                         GOTO 15                         
  604.    40 AL14 = ALPH1/4.0
  605.       AL24 = ALPH2/4.0
  606.       PT14 = PT(AL14,NDF)
  607.       PT24 = PT(AL24,NDF)
  608.       X0  = SQRT(C)                                                     
  609.       X1 = PT14*X0                                                      
  610.       X2 = PT24*X0                                                      
  611.    45 GM12 = GMA1/2.0
  612.       GM22 = GMA2/2.0
  613.       Z1 = PNOR(GM12)
  614.       Z2 = PNOR(GM22)
  615.       ALP12 = 1.0-ALPH1/2.0
  616.       ALP22 = 1.0-ALPH2/2.0
  617.       CHI1 = PCHI2(ALP12,NDF)
  618.       CHI2 = PCHI2(ALP22,NDF)
  619.       X1 = SQRT(VE)*(X1 + Z1*SQRT(NDF/CHI1))                            
  620.       X2 = SQRT(VE)*(X2 + Z2*SQRT(NDF/CHI2))                            
  621.                                         GOTO 60                         
  622.    50 AL12 = ALPH1/2.0
  623.       AL22 = ALPH2/2.0
  624.       NNDF = N-NDF
  625.       PF1 = PF(AL12,NNDF,NDF)
  626.       PF2 = PF(AL22,NNDF,NDF)
  627.       X1 = SQRT((N-NDF)*PF1*C)                                          
  628.       X2 = SQRT((N-NDF)*PF2*C)                                          
  629.                                         GOTO 45                         
  630. C                                                                       
  631.    60 YL1 = YES - X1                                                    
  632.       YL2 = YES - X2                                                    
  633.       YU1 = YES + X1                                                    
  634.       YU2 = YES + X2                                                    
  635. C                                                                       
  636.       RETURN                                                            
  637.       E N D                                                             
  638. C
  639.       SUBROUTINE   MINMAX(A,N,AMIN,AMAX)
  640. C
  641. C     ** Minimum and Maximum Value of Vector **
  642. C
  643. C     * REFERENCE, T.Haga & S.Hashinoto: ó╢▓╖╠▐▌╛╖ ─ ╝¡╛▓╠▐▌╠▐▌╛╖ú, P.20
  644. C
  645.       REAL*4       A(1)
  646. C
  647.       AMIN = A(1)
  648.       AMAX = A(1)
  649. C
  650.       DO 10 I=1,N
  651.          IF (A(I).LT.AMIN) AMIN = A(I)
  652.          IF (A(I).GT.AMAX) AMAX = A(I)
  653.    10 CONTINUE
  654. C
  655.       RETURN
  656.       END
  657. C
  658.       SUBROUTINE   DWRTIO(X,N,DWR)
  659. C
  660. C     ** Durbin-Watson Ratio **
  661. C
  662. C     * REFERENCE, T.Haga & S.Hashimoto: ó╢▓╖╠▐▌╛╖ ─ ╝¡╛▓╠▐▌╠▐▌╛╖ú, P.189
  663. C
  664.       REAL*4       X(1)
  665. C
  666.       AX = X(1)
  667.       SX = 0.0
  668.       DW = 0.0
  669. C
  670.       DO 10 I=2,N
  671.          FI = I
  672.          D = X(I)-AX
  673.          AX = AX+D/FI
  674.          SX = SX+D**2*(FI-1.0)/FI
  675.          DW = DW+(X(I)-X(I-1))**2
  676.    10 CONTINUE
  677. C
  678.       DWR = DW/SX
  679. C
  680.       RETURN
  681.       END
  682.