home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 18 REXX / 18-REXX.zip / mathapps.zip / SurfaceFit1.cmd < prev   
OS/2 REXX Batch file  |  1998-08-28  |  26KB  |  684 lines
  1. /* Typing "SurfaceFit1" without parameters will provide a help screen.       */
  2. /*                                                                           */
  3. /* Read a file of X1, X2, X3, ..., Y values and obtain a least squares       */
  4. /* estimate of a first degree surface fit to the data.                       */
  5. /*                                                                           */
  6. /* If this routine is called as an external routine a normal termination     */
  7. /* will return a value of 1.  A 0 is an error return.                        */
  8. /*                                                                           */
  9. /* The program makes use of the REXXLIB routines.  It is left to the user to */
  10. /* obtain and install this shareware package.  It is readily availble at     */
  11. /* Hobbes and other ftp sites, as well as from Quercus.                      */
  12. /*                                                                           */
  13. /* I wish to thank Mike Montenegro for his criticisms and assistance!        */
  14. /*                                                                           */
  15. /* This program may be copied and modified as desired provided appropriate   */
  16. /* credit is given.  Hopefully others will see fit to improve on it.         */
  17. /*                                                                           */
  18. /* Doug Rickman April 13, 1998 Marshall Space Flight Center, NASA.           */
  19. /*                                                                           */
  20. /* Programmer's notes:                                                       */
  21. /* I have done only one thing that is at all tricky, the rest is very ho-hum.*/
  22. /* I wanted the linear algebra routines to be generic subroutines, so they   */
  23. /* could be called and used in other programs.  To do this they must be able */
  24. /* to accept arrays of any name.  Since a subroutine in REXX can not return  */
  25. /* an array or even a compound variable the array must be created by the     */
  26. /* calling routine and the subroutine must work with it.  The Gauss-Jordan   */
  27. /* routine does its matrix operations in place!  In other words it writes    */
  28. /* the solution into the original, source matrices.  The two facts caused me */
  29. /* to set up the subroutines as follows -                                    */
  30. /*       The calling routine sets a variable called "elist" (for "edit list")*/
  31. /*    prior to calling the subroutine.  elist will contain the names of all  */
  32. /*    the matrices and other parameters needed by the subroutine.   The      */
  33. /*    subroutine is then called using only elist as a parameter.             */
  34. /*       In the subroutine elist is opened and the external variable names   */
  35. /*    are set to internal variable names.  The routine then uses the         */
  36. /*    "interpret" instruction to get the name of and work with the external  */
  37. /*    variable name.                                                         */
  38. /*                                                                           */
  39. /* This is decidedly a pain.  I could have also accomplished the objective   */
  40. /* using queues but I thought this might be quicker.  There is a decided     */
  41. /* performance penalty, but at this time it is not a real issue.             */
  42. /*                                                                           */
  43. /* I have written this using abundant comments and white space.  In fact     */
  44. /* approximately 1/3 of the total lines are "waste."  I hope this will help  */
  45. /* others who wish to edit the code.                                         */
  46. /*                                                                           */
  47.  
  48. signal on Halt
  49. signal on NotReady
  50.  
  51. Numeric Digits 9       /* This may be changed to increase/decrease precision */
  52.  
  53. if rxfuncquery('rexxlibregister') then do         /* this will start rexxlib */
  54.     call rxfuncadd 'rexxlibregister', 'rexxlib', 'rexxlibregister'  
  55.     call rexxlibregister
  56.     end
  57.  
  58. if rxfuncquery('sysloadfuncs') then do           /* this will start rexxutil */
  59.     CALL RxFuncAdd 'SysLoadFuncs', 'RexxUtil', 'SysLoadFuncs' 
  60.     CALL SysLoadFuncs
  61.     end
  62.  
  63. arg in print
  64. in    =strip(in)
  65. degree=strip(degree)
  66. print =strip(print)
  67.  
  68. if in='' | in='?' | in='-?' | in='/?' then call Help
  69. if dosisfile(in)<>1 then do
  70.    say 'The input file: ' in' is not a valid file.'
  71.     exit
  72.     end /* do */
  73.  
  74. /* --------------------------------------------------------------------------*/
  75. /* --- begin MAIN                                               -------------*/
  76. /* Read the data                                                             */
  77. /* Data are assumed to be x1, x2, x3, ..., y pairs, one per line, with or    */
  78. /* without delimiting commas.                                                */
  79. rc=stream(in,'c','open read')
  80. data=linein(in)
  81. NIndependents=words(data)-1
  82. Xs=''
  83. do i=1 to NIndependents
  84.    Xs=Xs 'x.'i'.n'
  85.    end i
  86. rc=stream(in,'C','seek =' 1)
  87.  
  88. n=0
  89. do while lines(in)\=0
  90.    n=n+1
  91.    data=linein(in)
  92.    data=translate(data,' ',',')
  93.    interpret 'parse var data 'Xs' y.n'
  94.    x.0.n=1 /* Used in the Sums subroutine */
  95.    end
  96. x.0=n
  97. rc=lineout(in)
  98.  
  99. /* Compute the individual sums */
  100. call Sums
  101.  
  102. /* Build the two matrices */
  103. do i=0 to NIndependents
  104.    ii=i+1
  105.    do j=0 to NIndependents
  106.       jj=j+1
  107.       A.ii.jj=SumX.i.j
  108.       end j
  109.    end i
  110. NRowsA=NIndependents+1
  111. NColsA=NIndependents+1
  112.  
  113. do i=0 to NIndependents
  114.    ii=i+1
  115.    b.ii.1=SumYX.i
  116.    end i
  117. NColsB=1
  118. NRowsB=NIndependents+1
  119.  
  120. /* Backup the original A array. */
  121. do i=1 to NRowsA
  122.    do j=1 to NRowsA
  123.       AA.i.j=A.i.j
  124.       end j
  125.    end i
  126.  
  127. /* now call the Gauss-Jordan routine */
  128. elist= 'A. NRowsA NColsA B. NRowsB NColsB'
  129. call GaussJordan elist
  130.  
  131. /* Critical quality check!!!!!                                               */
  132. /* Multiply original A. by A inverse to check for identity.                  */
  133. elist= 'A. NRowsA NRowsB AA. NRowsA NColsA C. NRowsC NColsC'
  134. call MatrixMultiply    elist
  135.  
  136. elist= 'C. NRowsC NColsC'
  137. rc=IndentityCheck()
  138. if rc=1 then do
  139.    /* say 'The result is close to the identity matrix.  All is well.'  */
  140.    end
  141. else do
  142.    say 
  143.    say 'Checking A x A(inverse):  Oh Oh!!'
  144.    say 
  145.    say 'The inverse of A. may not be numerically precise enough.  You should '
  146.    say 'examine the values in A. x A.inverse, in theory this should be the '
  147.    say 'identity matrix.  You may need to either increase the precision of '
  148.    say 'Numeric Digits (set at the beginning of the program) or increase the '
  149.    say 'amount of Numeric Fuzz set in the subroutine IndentityCheck:.  Of '
  150.    say 'course you might also consider putting in better data.  :-)   You can'
  151.    say 'also change the value of the variable "QualityChecks" internally in '
  152.    say 'the code to have the software dump additional information.'
  153.    say
  154.    end
  155.  
  156. /* Optional quality checks.  Set "QualityChecks=YES" to use.                 */
  157. QualityChecks='NO'
  158. if QualityChecks=YES then call QualityChecks
  159.  
  160. /* Print the solution matrix */ /* There are two formats available */
  161. say
  162. say 'Solution: '
  163.  
  164. /* This format will print the solution with ascending exponents */
  165. Solution=b.1.1
  166. do i=2 to NIndependents+1
  167.    Solution=Solution' + 'b.i.1'*x.'i-1
  168.    end i
  169. say 'y = 'Solution
  170.  
  171. /* Determine the standard error */
  172. call StandardError
  173.  
  174. /* To print the table or not to print! That is the question. */
  175.  
  176. select
  177.    when print='Y' then call ScreenPrint
  178.    when print='N' then nop
  179.    when validname(print) then do
  180.       if dosisfile(print) then do
  181.          say 'The file ' print' already exists, do you want it overwritten?'
  182.           say 'Enter a "y" for yes, "h" will give help,'
  183.          say 'any other response will abort processing.'
  184.          key=translate(sysgetkey())
  185.          say
  186.          select
  187.             when key='Y' then do
  188.                     rc=dosdel(print)
  189.                     if rc=1 then nop
  190.                     else do
  191.                        say 'The file 'print' could not be deleted.  Goodbye, Sweet Prince.'
  192.                         exit
  193.                         end /* else do */
  194.                    end /* when key='Y' then ... */
  195.  
  196.                 when key='H' then call Help
  197.                 otherwise exit
  198.                 end /* select */
  199.  
  200.             end /* if dosisfile(print) */
  201.       call FilePrint
  202.       end /* when validname(print) ... */
  203.    otherwise do
  204.       say
  205.       say 'Do you wish to see a table of x, y, y estimate, delta, delta^2?'
  206.       say 'Enter a "y" for yes, any other response will exit the program.'
  207.       key=translate(sysgetkey())
  208.       say
  209.       if key='Y' then call ScreenPrint
  210.       end
  211.    end /* select */
  212.  
  213. return 1
  214.  
  215. /* --- end MAIN                                                 -------------*/
  216. /* --------------------------------------------------------------------------*/
  217.  
  218. /* --------------------------------------------------------------------------*/
  219. /* --- begin subroutine - Help:                                 -------------*/
  220. Help:
  221. rc= charout(,'1b'x||'[31;7m'||'SurfaceFit1:'||'1b'x||'[0m')
  222. say ' Fit a 1st degree polynomial (a plane) to x1, x2, x3, ..., y data '
  223. say '       by the method of least squares.'
  224.  
  225. rc= charout(,'1b'x||'[33;1m'||'usage:'||'1b'x||'[0m')
  226. say ' SurfaceFit1 in print'
  227.  
  228. rc= charout(,'1b'x||'[33;1m'||'where:'||'1b'x||'[0m')
  229. say ' in     = the file holding x1, x2, x3, ..., y pairs, one set per line,'
  230. say '                comma or blank delimited.'
  231. say '       print  = if "y" a table of x, y, y estimate, delta, delta^2 is printed'
  232. say '                to the screen.'
  233. say '              = if "n" then the table will not be printed.'
  234. say '              = if a syntactically correct file name, print to this file.'
  235. say ''
  236.  
  237. rc= charout(,'1b'x||'[33;1m'||'Exam: '||'1b'x||'[0m')
  238. say ' "surfacefit1 Surface.data surface.out"  or "surfacefit1 Surface.data y"'
  239.  
  240. rc= charout(,'1b'x||'[33;1m'||'notes:'||'1b'x||'[0m')
  241. say ' If the parameter "print" is not given on the command line the user will'
  242. say 'be asked if the table is to be printed.  Output will go to the screen.'
  243. say 'If a file is specified only syntax is checked.  It is assumed the directory'
  244. say 'exists.  When printing to a file the full precision of the data are retained'
  245. say 'This can make for rather wide lines.  :-)  Precision is increased by changing'
  246. say 'the value NUMERIC DIGITS (a min. of 9 is assumed in the FilePrint subroutine).'
  247. say 'WARNINGS: There are no checks on user input.  For example, if the input '
  248. say 'file has a blank line SurfaceFit1 will probably go Boom.  There is a check of'
  249. say 'the inverse matrix computed (an essential step).  If quality is not very'
  250. say 'high the user is notified and suggestions made to address the problem.'
  251.  
  252. say 'Doug Rickman  April 20, 1998'
  253. exit
  254. return
  255.  
  256. /* --- end  subroutine - Help:                                  -------------*/
  257. /* --------------------------------------------------------------------------*/
  258.  
  259. /* --------------------------------------------------------------------------*/
  260. /* --- begin subroutine - Sums:                                 -------------*/
  261. /* Computes the sum of X, Y, X^2, Y^2, .....                                 */
  262. Sums:
  263.  
  264. /* Initiallize */
  265. do i=0 to NIndependents
  266.    do j=0 to NIndependents
  267.       SumX.i.j=0
  268.       end j
  269.    SumYX.i=0
  270.    end i
  271.  
  272. /* Sum of X and X1*X2 ... */
  273. do i=0 to NIndependents
  274.    do n=1 to x.0
  275.       do j=0 to NIndependents
  276.          SumX.i.j=SumX.i.j+x.i.n*x.j.n
  277.          end j
  278.       SumYX.i=SumYX.i+x.i.n*y.n
  279.       end n
  280.    end i
  281.  
  282. /* To list the sums, remove the comment marks. */
  283. /*
  284. do i=0 to NIndependents
  285.    do j=0 to NIndependents
  286.       say 'SumX.'i'.'j'='SumX.i.j
  287.       end j
  288.    say 'SumYX.'i'='SumYX.i
  289.    end i
  290. */
  291. return
  292. /* --- end subroutine   - Sums:                                 -------------*/
  293. /* --------------------------------------------------------------------------*/
  294.  
  295. /* --------------------------------------------------------------------------*/
  296. /* --- begin subroutine - GaussJordan:                          -------------*/
  297. /* Solve a square array using the Gauss-Jordan algorithm as outlined in      */
  298. /* numerical recipes.                                                        */
  299. /* Implemented by Doug Rickman March 13, 1998                                */
  300. /* elist holds:                                                              */
  301. /*      name of the first array,                                             */
  302. /*      the variable holding the number of rows in the first array,          */
  303. /*      the variable holding the number columns in the first array.          */
  304. /*      name of the second array,                                            */
  305. /*      the variable holding the number of rows in the second array,         */
  306. /*      the variable holding the number columns in the second array.         */
  307. /*                                                                           */
  308. /*      A. NRowsA NColsA B. NRowsB NColsB                                    */
  309. GaussJordan:
  310. procedure expose (elist)
  311.  
  312. parse var elist VArrayName1 VNRows1 VNCols1 VArrayName2 VNRows2 VNCols2 
  313.  
  314. NRows1=value(VNRows1)
  315. NCols1=value(VNCols1)
  316. ArrayName1=strip(VArrayName1,'t','.')
  317.  
  318. NRows2=value(VNRows2)
  319. NCols2=value(VNCols2)
  320. ArrayName2=strip(VArrayName2,'t','.')
  321.  
  322. N=NRows1     /* number of elements and rows */
  323. M=NCols2     /* right hand vactors is an array N by M */
  324.  
  325. NMax=50
  326. do j=1 to N
  327.    IPIV.j=0
  328.    end j
  329.  
  330. do i=1 to N
  331.    BIG=0
  332.    do j=1 to N
  333.       if IPIV.j <> 1 then do k=1 to N
  334.          interpret 'Temp='Arrayname1'.j.k'
  335.          if IPIV.k=0 & abs(Temp)>=BIG then do
  336.             interpret 'BIG=abs('Arrayname1'.j.k)'
  337.             irow=j
  338.             icol=k
  339.             end 
  340.          else if IPIV.k > 1 then do
  341.             say 'Singular matrix! Stop 1'
  342.             exit
  343.             end
  344.          end k
  345.       end j
  346.  
  347.    IPIV.icol=IPIV.icol+1
  348.  
  349.    if irow<>icol then do
  350.       do /*14*/ L=1 to N
  351.          interpret 'DUM='Arrayname1'.irow.L'
  352.          interpret Arrayname1'.irow.L='Arrayname1'.icol.L'
  353.          interpret Arrayname1'.icol.L=DUM'
  354.          end L
  355.       do L=1 to M
  356.          interpret 'DUM='Arrayname2'.irow.L'
  357.          interpret Arrayname2'.irow.L='Arrayname2'.icol.L'
  358.          interpret Arrayname2'.icol.L=DUM'
  359.          end L
  360.       end /* if irow<>icol then do ... */
  361.  
  362.    INDXR.i=irow
  363.    INDXC.i=icol
  364.    interpret 'Temp='Arrayname1'.icol.icol'
  365.    if Temp=0 then do
  366.       say 'Singular matrix! Stop 2.'
  367.       exit
  368.       end
  369.    
  370.    interpret 'PIVINV=1/'Arrayname1'.icol.icol'
  371.    interpret Arrayname1'.icol.icol=1'
  372.    do L=1 to N
  373.       interpret Arrayname1'.icol.L='Arrayname1'.icol.L*PIVINV'
  374.       end L
  375.    do L=1 to M
  376.       interpret Arrayname2'.icol.L='Arrayname2'.icol.L*PIVINV'
  377.       end L
  378.    do LL=1 to N
  379.       if LL \= icol then do
  380.          interpret 'DUM='Arrayname1'.LL.icol'
  381.          interpret Arrayname1'.LL.icol=0'
  382.          do L=1 to N
  383.             interpret Arrayname1'.LL.L='Arrayname1'.LL.L-'Arrayname1'.icol.L*DUM'
  384.             end L
  385.          do L=1 to M
  386.             interpret Arrayname2'.LL.L='Arrayname2'.LL.L-'Arrayname2'.icol.L*DUM'
  387.             end L
  388.          end /* if LL \= icol then do ... */      
  389.       end LL
  390.  
  391.    end i
  392.  
  393. /* Unscramble */
  394. do L=N to 1 by -1
  395.    if INDXR.L \= INDXC.L then do K=1 to N
  396.       INDXRL=INDXR.L
  397.       INDXCL=INDXC.L
  398.       interpret 'DUM='Arrayname1'.K.INDXRL'
  399.       interpret Arrayname1'.K.INDXRL='Arrayname1'.K.INDXCL'
  400.       interpret Arrayname1'.K.INDXCL=DUM'
  401.       end K
  402.    end L
  403.  
  404. return
  405. /* --- end subroutine   - GaussJordan:                          -------------*/
  406. /* --------------------------------------------------------------------------*/
  407.  
  408. /* --------------------------------------------------------------------------*/
  409. /* --- begin subroutine - MatrixMultiply:                       -------------*/
  410. /* Returns a 1 if successful, result is in C.                                */
  411. /* Returns a 2 if the matrices are not the correct size.                     */
  412. /* elist holds:                                                              */
  413. /*      name of the first array,                                             */
  414. /*      the variable holding the number of rows in the first array,          */
  415. /*      the variable holding the number columns in the first array.          */
  416. /*      name of the second array,                                            */
  417. /*      the variable holding the number of rows in the second array,         */
  418. /*      the variable holding the number columns in the second array.         */
  419. /*      name of the result array,                                            */
  420. /*      the variable holding the number of rows in the result array,         */
  421. /*      the variable holding the number columns in the result array.         */
  422. /*                                                                           */
  423. /*      A. NRowsA NColsA B. NRowsB NColsB C. NColsC NRowsC                   */
  424.  
  425. MatrixMultiply:
  426. procedure expose (elist)
  427.  
  428. parse var elist VArrayName1 VNRows1 VNCols1 ,
  429.                 VArrayName2 VNRows2 VNCols2 ,
  430.                 VArrayNameR VNRowsR VNColsR
  431.  
  432. NRows1=value(VNRows1)
  433. NCols1=value(VNCols1)
  434. ArrayName1=strip(VArrayName1,'t','.')
  435.  
  436. NRows2=value(VNRows2)
  437. NCols2=value(VNCols2)
  438. ArrayName2=strip(VArrayName2,'t','.')
  439.  
  440. ArrayNameR=strip(VArrayNameR,'t','.')
  441.  
  442. if NCols1=NRows2 then do
  443.    NRowsR=Nrows1
  444.    NColsR=NCols2
  445.    do i=1 to NRowsR
  446.       do j=1 to NColsR
  447.          interpret ArrayNameR'.i.j=0'
  448.          do k=1 to NCols1
  449.             interpret ArrayNameR'.i.j='ArraynameR'.i.j+('ArrayName1'.i.k*'ArrayName2'.k.j)'
  450.             end k
  451.          end j
  452.       end i
  453.    end /* if ... */
  454.  
  455.    interpret VNRowsR'=NRowsR'
  456.    interpret VNColsR'=NColsR'
  457.  
  458.    return 1
  459.    end
  460. else return 2
  461.  
  462. /* --- end subroutine   - MatrixMultiply:                       -------------*/
  463. /* --------------------------------------------------------------------------*/
  464.  
  465. /* --------------------------------------------------------------------------*/
  466. /* --- begin subroutine - IndentityCheck:                       -------------*/
  467. /* Check to see if matrix is an identity matrix (diagonal=1, other=0).       */
  468. /* 1 is returned for and identity matrix, a 0 otherwise.                     */
  469. /* Precision can be increased by increasing NUMERIC DIGITS for the Gauss-    */
  470. /* Jordan subroutine.  Also the tolerance for error can be increased by      */
  471. /* changing the NUMERIC FUZZ value in this routine.                          */
  472. /* elist holds:                                                              */
  473. /*          name of the array,                                               */
  474. /*          the variable holding the number of rows,                         */
  475. /*          the variable holding the number columns.                         */
  476. IndentityCheck: procedure expose (elist)
  477.  
  478. parse var elist ArrayName1 VNRows1 VNCols1
  479.  
  480. NRows=value(VNRows1)
  481. NCols=value(VNCols1)
  482. Arrayname=strip(Arrayname1,'t','.')
  483.  
  484. n=digits()
  485. numeric fuzz n-3
  486. do i=1 to NRows
  487.    do j=1 to NCols
  488.       interpret 'test=1+'Arrayname'.i.j'
  489.       if i=j & test=2 then iterate
  490.       if test=1 then iterate
  491.       else do
  492.          numeric fuzz 0
  493.          return 0
  494.          end 
  495.       end j
  496.    end i
  497. numeric fuzz 0
  498. return 1
  499. /* --- end subroutine   - IndentityCheck:                       -------------*/
  500. /* --------------------------------------------------------------------------*/
  501.  
  502. /* --------------------------------------------------------------------------*/
  503. /* --- begin subroutine - ShowMatrix:                           -------------*/
  504. /* elist holds:                                                              */
  505. /*          name of the array,                                               */
  506. /*          the variable holding the number of rows,                         */
  507. /*          the variable holding the number columns.                         */
  508.  
  509. ShowMatrix: procedure expose (elist)
  510.  
  511. parse var elist ArrayName1 VNRows1 VNCols1
  512.  
  513. NRows=value(VNRows1)
  514. NCols=value(VNCols1)
  515. Arrayname=strip(Arrayname1,'t','.')
  516. /* say 'Array: 'ArrayName */
  517. do i=1 to NRows
  518.    row=''
  519.    do j=1 to NCols
  520.       interpret 'row=row 'ArrayName'.i.j'
  521.       end j
  522.    say row
  523.    end i
  524. return
  525.  
  526. /* --- end subroutine   - ShowMatrix:                           -------------*/
  527. /* --------------------------------------------------------------------------*/
  528.  
  529. /* --------------------------------------------------------------------------*/
  530. /* --- begin subroutine -                                       -------------*/
  531. QualityChecks:
  532.  
  533. say
  534. /* Print the inverse of A. */
  535. say 'Inverse of A is '
  536. elist= 'A. NRowsA NColsA'
  537. call ShowMatrix elist
  538.  
  539. say
  540. /* Multiply A inverse times the solution vector.  The result should equal the*/
  541. /* original B. matrix.                                                       */
  542. elist= 'AA. NRowsA NRowsB B. NRowsB NColsB C. NRowsC NColsC'
  543. say 'Inverse of A times the solution B: (should equal the original matrix B.)'
  544. call MatrixMultiply    elist
  545. elist= 'C. NRowsC NColsC'
  546. call ShowMatrix elist
  547.  
  548. say
  549. /* Multiply original by inverse to check for identity.                       */
  550. elist= 'A. NRowsA NRowsB AA. NRowsA NColsA C. NRowsC NColsC'
  551. say 'Original A times the inverse of A: (should equal identity matrix)'
  552. call MatrixMultiply    elist
  553. elist= 'C. NRowsC NColsC'
  554. call ShowMatrix elist
  555.  
  556. return
  557.  
  558. /* --- end subroutine   -                                       -------------*/
  559. /* --------------------------------------------------------------------------*/
  560.  
  561. /* --------------------------------------------------------------------------*/
  562. /* --- begin subroutine - StandardError:                        -------------*/
  563. StandardError:
  564.  
  565. SE=0
  566. do n=1 to x.0
  567.    yHat=b.1.1
  568.    do i=2 to NIndependents+1
  569.       ii=i-1
  570.       yHat=yHat + b.i.1*x.ii.n
  571.       end i
  572.    SE=SE+(y.n-yHat)*(y.n-yHat)
  573.    end n
  574. SE=sqrt(SE/x.0)
  575.  
  576. say
  577. say 'Standard Error of Estimate = 'SE
  578.  
  579. return
  580. /* --- end  subroutine  - StandardError:                        -------------*/
  581. /* --------------------------------------------------------------------------*/
  582.  
  583. /* --------------------------------------------------------------------------*/
  584. /* --- begin subroutine - ScreenPrint:                           -------------*/
  585. ScreenPrint:
  586. Header1=''
  587. Header2=''
  588. Table=''
  589. do i=1 to NIndependents
  590.    Header1=Header1'    x.'i' '
  591.    Header2=Header2'_______ '
  592.    Table=  Table'left(x.'i'.n,7) '
  593.    end i
  594.  
  595. say
  596. say Header1 '     y y estimate    delta    delta^2'
  597. say Header2 '______ __________ __________ __________'
  598.  
  599. do n=1 to x.0
  600.    yHat=b.1.1
  601.    do i= 1 to NIndependents
  602.       ii=i+1
  603.       yHat=yHat + b.ii.1*x.i.n
  604.       end i
  605.   
  606.    interpret 'say ' Table ,
  607.                 ||'right(y.n,7)',
  608.                 ||' right(yHat,10)',
  609.                 ||' right(format(y.n-yHat,,3,3,3),10)',
  610.                 ||' right(format(pow(y.n-yHat,2),,3,3,3),10)'
  611.    end n
  612.  
  613. return
  614. /* --- end  subroutine  - ScreenPrint:                          -------------*/
  615. /* --------------------------------------------------------------------------*/
  616.  
  617. /* --------------------------------------------------------------------------*/
  618. /* --- begin subroutine - FilePrint:                            -------------*/
  619. FilePrint:
  620.  
  621. rc=lineout(print,'Input file: 'in)
  622. rc=lineout(print,'Date:       'date())
  623. rc=lineout(print,'Time:       'time())
  624. Solution=b.1.1
  625. do i=2 to NIndependents+1
  626.    Solution=Solution' + 'b.i.1'*x.'i-1
  627.    end i
  628. rc=lineout(print,'y = 'Solution)
  629. rc=lineout(print,'Standard Error of Estimate: 'SE)
  630. rc=lineout(print,'')
  631.  
  632. Precision=digits()+3
  633. Header1=''
  634. Header2=''
  635. Table=''
  636.  
  637. Header1=right('n',Precision,' ')
  638. Header2=right('',Precision,'_')
  639.  
  640. do i=1 to NIndependents
  641.  
  642.    Header1=Header1 right('x.'i,Precision,' ')
  643.    Header2=Header2 right('',Precision,'_')
  644.    Table=  Table 'right(x.'i'.n,Precision) '
  645.    end i
  646.  
  647. say
  648. rc=lineout(print,Header1 right('y',Precision) right('y estimate',Precision),
  649.            right('delta',Precision) right('delta^2',Precision))
  650. rc=lineout(print,Header2 right('',Precision,'_') right('',Precision,'_') ,
  651.            right('',Precision,'_') right('',Precision,'_'))
  652.  
  653. do n=1 to x.0
  654.    yHat=b.1.1
  655.    do i= 1 to NIndependents
  656.       ii=i+1
  657.       yHat=yHat + b.ii.1*x.i.n
  658.       end i
  659.   
  660.    interpret 'rc=lineout(print,right(n,Precision)' Table ,
  661.                 ||' right(y.n,Precision)',
  662.                 ||' right(yHat,Precision)',
  663.                 ||' right(format(y.n-yHat,,3,3,3),Precision)',
  664.                 ||' right(format(pow(y.n-yHat,2),,3,3,3),Precision))'
  665.    end n
  666.  
  667. return
  668. /* --- end  subroutine  - FilePrint:                            -------------*/
  669. /* --------------------------------------------------------------------------*/
  670.  
  671. /* --------------------------------------------------------------------------*/
  672. /* --- begin subroutine - Halt:                                 -------------*/
  673. Halt:
  674. say 'This is a graceful exit from a Cntl-C'
  675. return 0
  676. /* --- end  subroutine - Halt:                                  -------------*/
  677. /* --------------------------------------------------------------------------*/
  678. /* --- begin subroutine - NotReady:                             -------------*/
  679. NotReady:
  680. say 'It would seem that you are pointing at non-existant data.  Oops.  Bye!'
  681. return 0
  682. /* --- end  subroutine - NotReady:                              -------------*/
  683. /* --------------------------------------------------------------------------*/
  684.