home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / sigm / vol269 / regrs.doc < prev    next >
Encoding:
INI File  |  1986-05-22  |  13.3 KB  |  281 lines
  1. [REGRS.DOC of JUGPDS Vol.12]
  2.  
  3.       ** MULTIPLE LINEAR REGRESSION ANALYSIS WITH INTERVAL ESTIMATION **
  4.                  (╛▌╣▓ ╝▐¡│╢▓╖ ╠▐▌╛╖ ─ ╕╢▌╜▓├▓)
  5.  
  6.     Yoshio MONMA (JUG-CP/MS No.43)            85-03-21    
  7.     1514-63 Takata-cho, Khohoku-ku
  8.     Yokohama 223, Japan
  9.         Phone: 03-719-2271 (office), 045-543-5232 (home)
  10.  
  11. 1. INTRODUCTION
  12.     The multiple regression analysis is one of statistical methods 
  13. commonly used in various applications.    Although this technique is well-
  14. developed for the mainframe machines from early stage of computer usage,
  15. here is a version for CP/M on the micro.  The program has been written 
  16. in Microsoft's FORTRAN-80 (Version 3.44) under CP/M-80 V.2.2.
  17.  
  18.     Due to the restriction of memory space available in CP/M, I had
  19. to split the program into two parts: REGRS1.FOR(REG1.COM) and 
  20. REGRS2.FOR(REG2.COM).  The main job of input data, the regression 
  21. procedure, the analysis of variance, and the interval estimation are
  22. processed by REGRS1, while REGRS2 plots the scatter diagrams of the
  23. residuals to the LST: device (printer).
  24.  
  25. 2. THEORY
  26.     Since there are many good text books in this field [1-5],
  27. I won't repeat the details of mathematical manipulation which will
  28. require several pages to be complete.  Instead only minimal informa-
  29. tion needed to use the present program is given in the followings.
  30.  
  31.     The present version is largely based on the programs shown in
  32. the book titled ó╢▓╖╠▐▌╛╖─ ╝¡╛▓╠▐▌╠▐▌╛╖ú ("Regression Analysis and 
  33. Principal Component Analysis") by T. Haga and S. Hashimoto [1].  Those
  34. who want the general idea and the theoretical background should get 
  35. this book.  What I have added here to the original version by them,
  36. which was written in the conventional Fortran language for the 
  37. machines of minicomputer or above, is the implementation of statis-
  38. tical interval estimation (─│╣▓├╖ ╕╢▌ ╜▓├▓).  Also I owe it to 
  39. Professor T. Sawa through his exellent text book [2].
  40.  
  41.     The data matrix looks like as follows:
  42.                                                                        
  43.                1      2          J         JP    JPY  JPY1             
  44.               X1     X2   ...   XJ   ...   XP   XP+1  XP+2             
  45.          1  X(1,1) X(1,2) ... X(1,J) ... X(1,P) Y(1)   1               
  46.          2  X(2,1) X(2,2) ... X(2,J) ... X(2,P) Y(2)   1               
  47.          .    ...    ...  ...   ...  ...   ...   ..    .               
  48.          .    ...    ...  ...   ...  ...   ...   ..    .               
  49.          I  X(I,1) X(I,2) ... X(I,J) ... X(I,P) Y(I)   1               
  50.          .    ...    ...  ...   ...  ...   ...   ..    .               
  51.          .    ...    ...  ...   ...  ...   ...   ..    .               
  52.         IN  X(N,1) X(N,2) ... X(N,J) ... X(N,P) Y(N)   1               
  53.                                                                        
  54. where JP is the number of independent variables, and IN is the
  55. data (sample) size.
  56.  
  57.     For I=1,IN data points the multiple regression model is:
  58.                                                                       
  59.              Y(I) = B0 + B1*X(I,1) + B2*X(I,2) + ... + BP*X(I,P),      
  60.                                                                        
  61. where B0,B1,B2,...BP are the partial regression coefficients.     
  62.                                                                        
  63. 3. EXPLANATION OF THE PROGRAM
  64.    3.1    Nomenclature 
  65.  
  66.      IN           No. of Samples (n)
  67.      JPY       No. of Variables (p+1) 
  68.      MPL       Degree of Polynomial Regression (When JPY=2)
  69.      INP0      File reference no. for data input device (FORT0X.DAT)
  70.      INP1      File reference no. for X(I,J), usually INP1=INP0
  71.                INP0 and INP1 must be between 6 and 9 (FORT06.DAT - 
  72.                FORT09.DAT.
  73.      KDTA      When KDTA=1 no out put of the raw data
  74.      ISWP      If ISWP=1, the computation process is output
  75.      IRES      When IRES=1, no output of residual table and graph
  76.      JGRH      If JGRH=1, graphs of residuals vs. X(J) are output.
  77.                                                                        
  78.      INTV       Index of Interval Estimation:
  79.           = 0   Average (Mean Estimate) only,
  80.           = 1   Avg. + Confidence Intervals (CI),                      
  81.          ╝▌╫▓-╕╢▌ = ─╕├▓ ╔ ╛┬╥▓-═▌╜│┴ ╞ └▓╡│╜┘ ╝▐¡│┐▐╕-═▌╜│ (Y) ╔      
  82.                     ╖└▓┴ ╢▐ ICFL(%) ╔ ╝▌╫▓─▐ ├▐ ╠╕╧┌┘ ╪«│▓╖.           
  83.           = 2   Avg. + Simultaneous Conf. Intervals (SCI),             
  84.          ─▐│╝▐-╝▌╫▓-╕╢▌ = 1 ╕╨ ╔ ╦«│╬▌-╢▌┐╕┴ DATA ╢╫ 1┬ ╔ ╢▓╖╝╖ ª      
  85.               ╜▓├▓ ╝, ╠╕╜│║ ╔ ║─┼┘ ╛┬╥▓-═▌╜│ ╞ └▓╡│╜┘ Y ╔ ╖└▓┴ ╢▐      
  86.               ╝├▓╜┘ ╝▌╫▓─▐ ├▐ ╠╕╧┌┘ ╪«│▓╖.                             
  87.           = 3   Avg. + Prediction Intervals (PI),                      
  88.          ╓┐╕-╕╢▌ = ─╕├▓ ╔ ╛┬╥▓-═▌╜│ ╞ └▓╡│╜┘ Y ª 1╢▓ └▐╣ ╓┐╕ ╜┘.       
  89.           = 4   Avg. + Tolerance Intervals (TI),                       
  90.          ╖«╓│-╕╢▌ = ─╕├▓ ╔ ╛┬╥▓-═▌╜│┴ ╞ └▓╡│╜┘ ╓┐╕ ª ╕╪╢┤╝             
  91.                   ╡║┼│─╖, ┐║ ╞ ╠╕╧┌┘ ╢╕╪┬ ╢▐ ╜╕┼╕─╙ IPRB(%) ├▐▒┘ ─▓│║─ 
  92.                   ª ╝▌╫▓─▐ ICFL(%) ª╙»├ ╝¡┴«│ ├▐╖┘ ╩▌▓.                
  93.           = 5   Avg. + Simultaneous Tolerance Intervals (STI)          
  94.          ─▐│╝▐-╖«╓│-╕╢▌ = 1 ╕╨ ╔ ╦«│╬▌-╢▌┐╕┴ DATA ╞ ╙─┬▐╖, ╜▓├▓╝└      
  95.                   ╢▓╖╝╖ ª ╙┴▓├, ▓╕┬╢ ╔ ║─┼┘ ╛┬╥▓-═▌╜│┴ ╞ └▓╝├  
  96.                   ╝▌╫▓─▐ICFL(%) ╔ ╓┐╕ ª ┼▌╢▓ ╙ ╡║┼»└─╖, ┐╔     
  97.                   '╛▓║│╪┬' ╢▐ ▒╫╢╝▐╥ ╝├▓╗┌└ IPRB(%) ╞ ╦─╝▓╓│┼ ╕╢▌. 
  98.                                                                        
  99.      ICFL      Confidence level (%)                                   
  100.      IPRB      Probability level (%) for TI and STI                   
  101.      NPRD      No. of predictions to be made                          
  102.      CMNT       Any comments (18A4 = 72 characters)                    
  103.      SAMP      Sample/lot no. (A8)                                    
  104.      VARS      Description of variables (18A4)                        
  105.      FMT       Variable format (18A4)
  106.      X(I,J)    Data matrix                                            
  107.      XX(I,J)   Specified values of predictor variables to make
  108.            the estimation
  109.      ITRS(J)   Transformation code for variable X(I,J):               
  110.              = 0   No transformation
  111.              = 1   Linear: (X --> A*X+B)
  112.              = 2   Inversion: (X --> 1/X)
  113.              = 3   X --> LOG10(X)
  114.              = 4   X --> LOG(X)
  115.              = 5   X --> SQRT(X)
  116.              = 6   X --> X**2
  117.              = 7   X --> X**P
  118.              = 8   X --> 10.0**X
  119.              = 9   X --> EXP(X)
  120.              =10   Logistic: (X --> LOG(X/1.0-X))
  121.                                                                       
  122.      S(JP,JP)   Sum of squares
  123.      R(JP,JP)   Correlation matrix
  124.      JP         No. of independent variables
  125.      NDF        Degree of freedom (n-p-1 = IN-JP-1)                    
  126.      AVEX(J)    Average of X(J)                                        
  127.      SIGX(J)    Sigma of X(J)                                          
  128.      YEST(I)    Estimate of Y(I)                                       
  129.      RES(I)     Residuals                                               
  130.      TOL(J)     Tolerance to remove variables with high correlation
  131.  
  132.    3.2 Subprograms
  133.         Subprograms needed are:
  134.      BIGCHR, SSSP, ANOVA, ESTIM, MINMAX, DWRTIO  [REG1.FOR: This file]  
  135.      TRNSV0, MATPRI, VECPRI, CORREL, SWEEP1 [MATSUB.FOR]
  136.      PF, PT, PNOR, PCHI2 [PRBSUB.FOR]
  137.  
  138. BIGCHR    Print cndensed and enlarged string on RP-80
  139. TRNSV0    Transformation of a vector
  140. MATPRI    Print-out of a matrix  
  141. MINMAX    Minimum and maximum of a vector (one-dimensional array)
  142. VECPRI    Print out of a vector
  143. SSSP    Read the data matrix from FORT10.DAT and computes sum of squares
  144.     and sum of products
  145. ANOVA    Analysis of variance
  146. ESTIM    Estimation by the regression and residuals
  147. INTRVL    Interval estimation
  148. DWRTIO    Durbin-Watson ratio
  149. PT    Percentage point of t-distribution
  150. PF    Percentage point of F-distribution
  151. PNOR    Percentage point of the normal distribution
  152. PCHI2    Percentage point of Chi-squared distribution
  153.  
  154.    3.3 Structure of the program    
  155.     Here is an outline of the subroutine hierarchy for REGRS1.
  156.  
  157. REGRS1 (Main)
  158.     BIGCHR
  159.     TRNSV0
  160.     MATPRI
  161.     MINMAX
  162.     VECPRI
  163.     SSSP
  164.     ANOVA
  165.         PF
  166.         PT
  167.     ESTIM
  168.         PT
  169.         INTRVL
  170.             PF
  171.             PT
  172.             PNOR
  173.             PCHI2
  174.     DWRTIO
  175.  
  176.     As for REGRS2 it looks like:
  177.  
  178. REGRS2 (Main)
  179.     BIGCHR
  180.     MINMAX
  181.     GRAPH1
  182.  
  183.     The data matrix is read from FORT0X.DAT(INP1) and is saved to 
  184. FORT10.DAT.  The subroutine SSSP reads the data matrix from FORT10.DAT
  185. to make successive calculation of sums of squares and products with
  186. high accuracy.  The regression coefficients, residuals and interval
  187. estimation are also written to FORT10.DAT so that the data to plot
  188. the scatter diagrams of residuals is read by REGRS2.
  189.  
  190. 4. RESTRICTIONS AND HARDWARE DEPENDENCY
  191.     This program has been designed for Microsoft's FORTRAN-80 and 
  192. EPSON's RP-80 printer, therefore you must modify some statements 
  193. if you are going to run this program under another system.  I am going
  194. to implement this program for MS-Fortran.  It will be relatively easy 
  195. job, I think.  Here are some points:
  196.  
  197.         1) The data file must be PIPed before the execution because 
  198. this program uses the convention of defaulted file reference no. in 
  199. FORTRAN-80 from 6 (FORT06.DAT) to 9 (FORT09.DAT).  Examine the batch
  200. (*.SUB) files how to do this.  Don't use TAB code to make the data
  201. file I had experienced some difficult time to recognize the problem.
  202. The file FORT10.DAT is used to work and to save the data matrix and 
  203. results.
  204.  
  205.     2) Max. data size = 80(100) and  max. no. of independent 
  206. variables = 12(15) are dependent on the TPA size; 56k(64k).
  207.  
  208.     3) The printer routines BIGCHR and GRAPH are hardware dependent.
  209. You must be careful to check control code of your printer.  Both are
  210. not essential routines.
  211.  
  212.     I have checked this program with NEC PC-8001 Mk-II and 
  213. Kohjinsha's MICRO DECISION.  The latter machine was advantageous over
  214. the former because of larger TPA size (see above) and faster disk I/O.
  215.  
  216. 5. THE DATA FILE
  217.  ====================================================================  
  218.      ** List of Input Data (Card Image File) **
  219. ---------------------------------------------------------------------
  220.    FMT          Data (Variables)                        No. of cards  
  221. ---------------------------------------------------------------------
  222. (14I5)       IN,JPY,MPL,INP1,KDTA,ISWP,IRES,JGRH,                 
  223.               INTV,ICFL1,IPRB1,ICFL2,IPRB2,NPRD                       
  224.  (18A4,A8)    CMNT,SAMP                                        1       
  225.  (18A4)       VARS                                            JPY      
  226.  (18A4)       FMT                                              1       
  227.  FMT         ((X(I,J),J=1,JP),I=1,IN)                      variable    
  228.  ------------------------------------------------------------------       
  229.  [If NPRD<>0, read specified values of predictor var. ]               
  230.  (18A4)       FMT                                              1       
  231.  (FMT)        ((XX(I,J),J=1,JP),I=1,NPRD)                  variable    
  232.  -------------------------------------------------------------------      
  233.  (12I5)       (ITRS(J),J=1,JPY)                                1
  234.  --------------------------------------------------------------------        
  235.  [If ITRS(J)=1, read (A,B) for X = A*X+B]      
  236.  (2F10.3)          A,B                                         1       
  237.  [If ITRS(J)=7, read P for X = X**P]
  238.  (F10.2)           P                                           1
  239.  ====================================================================  
  240.  
  241.     Note that if 'TRNSV0' needs the parameters they must be supplied
  242. in the data file FORT0X.DAT.
  243.  
  244. 6. INSTALLATION
  245.    6.1 If you don't have FORTRAN-80 but available RP-80
  246.     Read carefully this document and prepare the data file in 
  247. such a way described in the previous section.  Run it just as shown
  248. in REGTST1.SUB and REGTST2.SUB.
  249.  
  250.    6.2 If you have both FORTRAN-80 and RP-80
  251.     You may recompile the source program with F80 and link with
  252. L80 to produce the COM files.  I have provided the complete command
  253. sequence to do this in REGRS.SUB.  Prior to issue the SUBMIT command
  254. you must copy the batch file to Drive A: (I think you know why) via:
  255.     A>PIP A:=B:*.SUB
  256. Note that you have your Fortran system (F80.COM and L80.COM) and the
  257. batch programs (SUBMIT.COM and XSUB.COM) on Drive A:, while the source
  258. code (*.FOR) must be in Drive B:.
  259.  
  260.     Then simply issue the batch command:
  261.     A>SUBMIT REGRS
  262. And if it seems to give satisfactory results, then try the other test 
  263. programs:
  264.     A>SUBMIT REGTST1
  265.     A>SUBMIT REGTST2 
  266.  
  267.    6.3 If you wish to modify or improve the program
  268.     Place your Fortran system and editor on Drive A:, *.FOR and 
  269. *.DAT files in Drive B:, respectively.  Do the cyclic loop of edit, 
  270. compile, link, and run until you are satisfied.
  271.  
  272.   
  273. REFERENCES                                                     
  274. [1] ╩╢▐Ñ╩╝╙─ (Haga and Hashimoto): ó╢▓╖╠▐▌╛╖ ─ ╝¡╛▓╠▐▌╠▐▌╛╖ú (REGRESSION
  275.     ANALYSIS AND PRINCIPAL COMPONENT ANALYSIS), ╞»╢╖▐┌▌ (JUSE) (1980)
  276. [2] ╗▄ └╢╨┬ (T. Sawa): ó╢▓╖╠▐▌╛╖ú (REGRESSION ANALYSIS), ▒╗╕╫╝«├▌ (1979)
  277. [3] ╡╕╔ ╬╢ (Okuno & others): ó└═▌╪«│╢▓╛╖╬│ú (MULTIVARIATE ANALYSIS),
  278.     Vol.1(1971), Vol.2(1976), ╞»╢╖▐┌▌ (JUSE).
  279. [4] M.G. Kendall: MULTIVARIATE ANALYSIS, Griffin, (1975).
  280. [5] 
  281.