Scien-Subroutines, Volume I (Bye/McGraw-Hill, Peterborough, NH, 1981).
  
 
345678
910qd.1213141516V718a920ai22232425262728293031323334353637383940al4243444646474849
APPENDIX 2 635Oo
SSTORAGE: 2PROGRAM LSTSOC---+ Least Squares Regression Analysis for 2 variables
CHARACTER*1 ANSDIMENSION X(100),¥¢100)
DATA MAX/100/OPEN(6, FILE=" PRN")
C---+ Read data in sets of 2 variablesTelWRITEC*, 7000)7000 FORMAT(' Enter data in sets of 2 (Y then X).'/' If second variable is omitted, sequential numbers’* will be used for X.'/" Extra <return> ends input. *//* Enter below:'/*OYYYY YY. XXXXXX")10 READ(*,5000)¥(1) ,XCI)5000, FORMAT(EN,F6.2,F7.2)IF(Y¥(1).£0.0)G0TO 30TF(X(1).€0.0)XCQ)=1IF (1.€0.MAX) THEN 'WRITE(*, 7001 MAX |7001, FORMAT(' Data input reached max of',15,'. Processing begins.")*  g0To 30ENDIF1 =I+1GOTO 10
seat
30. N= 1-1WRITE(*, 7002)7002 FORMAT(" Print data (Y/N)?"\)READ(*,5001)ANS5001 FORMAT(AL)IF(ANS.NE. 'Y*)GOTO 50WRITE(6, 6000)6000 FORMAT('1Least-Squares Regression Analysis’ ,20X,‘Date Input'\\+ " ¥ (dep) xX Cind)"\)90 40 J = 1,NWRITE(6 ,6001)¥(J), XC)6001 FORMAT(2F8.2)40 CONTINUE
50 CALL RAZV(Y,X,N) -
CALL EXIT ~END -
STITLE: 'RAZV: Regression for 2 Variables
SSUBTITLE: ‘Copypright 1986 PJ Kaufman’
SSTORAGE:SUBROUTINE RA2V(Y,X,N)
 
Regression Analysis for 2 VariablesWith Linear, Power and Log Transformations
Y = Dependent variable nee- X = Independent variableN= Number of data entries
  
636
1213
14
15
16
17
18
19
20
21 C----22
23
24
25
26
27
28 600029
30
31
32
33
34
35
36
37
38
39
40 60044
42
43
44
45
46
47
48 600149
50
51
52
53
54
55
56 600257
$3
59 700060
61 500062
63
64 700165 +86
67 $00163
69
70 600971
72 C-+--
CHARACTER*1 ANSCHARACTER*2 SIGN, SIGNCCHARACTER*9 YNAME, XNAME
DIMENSION X(100),¥(100),ZAC4),2B(4) ,2R(4)DATA MAX/100/
Perform linear regression
ITYPE=1
CALL BSTFIT(Y,X,N,ZA(1),28(1),Z2R(1),SD,1TYPE)SIGN="
IF(ZB(1).GE.0)SIGN=" +*
WRITE(6,6000)ZA(1), SIGN, Z8(1),2R(1),S0FORMAT(/54X, "2 St dev'//+ * LINEAR Y=", F8.3,A2,F8.3,'*X',18X,F7.1,F9.3)
CALL LSTSO2(¥,X;N,ZA(4),2B(4),C,S4)
SIGN ="Tf(ZB(4).GE.0)SIGN =" +°SIGNC ="1F(C.GE.O0)SIGNC’=" +"ITYPE = 4
WRITE(6 6004 )ZA(4) , SIGN, 28(4) ,SIGNC,C,S4FORMAT(/* CURVE Y=", F8.3,A2,F8.3,'*X’,A2,F8.3,1*X2",22X,F9.3)
ITYPE=2
CALL BSTFITCY,X,N,ZA(2),28(2),ZR(2), SD, ITYPE)SIGN=" '
TF(ZB(2).GE.0)SIGN=* +"
WRITE(6,6001)ZA(2), SIGN, 28(2),2R(2), SDFORMAT(/* LOG — Ln(¥)=",F8.3,A2,FB.3,‘*Ln(X)*,10X,F7.1,F9.3)
ITYPE=3
CALL BSTFITCY,X,N,ZA(3) ,2B(3) ,2R(3),SD,ITYPE)SIGN=" *
IF(ZB(3).GE.0)SIGN=' +°
WRITE( 6, 6002)ZA(3) , SIGN, ZB(3) ,ZR(3),SDFORMAT(/* EXPON Ln(Y¥)=",F10.3,A2,F10.3, '*X",10X,F7.1,F9=3)
uRITEC*, 7000)
FORMAT( Sample results (Y/N)?'\)
READ(* ,5000 ANS -
FORMAT(A1)
TFCANS. NE. "N" THEN -WRITEC*, 7001) “FORMAT(/' Enter start, end and increments below'/
"ssssss eeeeee iiiiti')READ(*,S001)P, END, STEPSFORMAT(BN,F6.2,2F7.2)
WRITE(6, 6009)FORMAT(//12X, 'X LINEAR CURVI 106 EXPON'/)
Calculate all four fits at the same time
 
 
 
APPENDIX 2 -73 100.0 YL = ZAC1) + ZBCL)*P74 Y2 = ZACA) + ZBC4)*P + CPAP75 ¥3 = ZA(2) + ZB(2Z)*ALOG(P)76 Y3 = EXP(Y3)7 Y4 = 2A(3) + 2B(3)*?78 Y4 = EXP(Y4)7980 WRITE(6,6010)P,¥1,¥2,¥3,¥481 6010 — FORMAT(SX,5F8.2)8283 1F(P.LT.END)THEN84 P= Pp + STEPS85 GOTO 10086 ENOIF8788 ENDIFag90 WRITEC*, 7002)91 7002 FORMAT(/' Plot results (¥/N)?"\)92 READ(*,5000) ANS93 TF(ANS.NE.'N*)THEN94 WRITE(*, 7003)95 7003 FORMAT(' Name for x-scale (9 chars)?")96 ' READ(*, 5002) XNAME97 5002  — FORMAT(A)98 WRITE(*,7004)99 7004 FORMAT(’ Name for y-scale (9 chars)?")100 READ(*, 5003) NAME301 $003 FORMAT(AS)102 120 © WRITE(*,7005)103. 7005  FORMAT(/" Which plot?'/’ 1: Linear'/* 2: Logarithmic'/104 + 3: Exponential’/* 4: Curvilinear'/105 + * Which number?"\)106 READ(*, 5004) 1107 5004 — FORMAT(BN, 14)108 Te =C109 IF(I.NE.4)TC = 0110 CALL PLOTRA(N;X,Y,ZACI),ZB(1),TC,1,KNAME, YNAME)ll. WRITE(*,7006)112 7006  FORMAT(/* Another plot (Y/N)?'\)113 READ(*,5000) ANS114 IFCANS.NE. ‘N")GOTO 120115 ENDIF116417 RETURN118 END -1 STITLE: 'LSTSQ2: 2nd Order Least Squeres’2 SSUBTITLE: ‘P.d. Kaufman’3 SSTORAGE:2 - _4 SUBROUTINE LSTSO2Z(¥,%,N,A,8,C,5)56 Second Order (Curvitinear) Regression for 2 variables78 Dependent variable9 Independent variable
 
noarzx<
Number of data entriesConstent term
Coefficient of X -Coefficient of xX**2
Standard deviation of residuals
637 
 
 
 
 
     
 
S38 APPENDIX? APPENDIX 2 : $39ae15 :i ---+ Transform data from log form to normal form based upan thej i$ DIMENSION X(2),¥C2) ; é mem type of regression fit (transmitted via ITYPE)18 DOUBLE PRECISION SXX,SKY,SYY,SXX2,SK2K2,SYX2 8 GoTo (110,120,130), 1TYPE: 20 C---- Find the mean of x and y q Ce--- Lineari 22 are 3 12.110 TY#A+B*x| . 23 DO 10 t= 1,N 3 GoTo 140J 24 AX = AX + X(T)15 C---- LogSe » AX a “rp 16 120 TY=A+B*ALOG(X}27 AY = AN 17 TYMEXP(TY)i 18 GoTo 1404 28 19{= 29 C---- Initialize for summations 20 C---- Exponentiai 1 sky = 0 21 130 TY=A+B¥xt : x svt 0, 2 TY = EXP(TY)t 34 Sx2x2 = 0 bs 140 REGVAL=TY38 SYK2 = 0 26 RETURN37 C---+ Sums 27 END38 00 20 1 = 1,N ' 11939 XT = x1) : 120 SUBROUTINE BSTFIT(Y,X,N,A,8,R, SD, 1TYPE}40 yD = YC) 12141 XM = XE - Ax : 122 € - Basic linear regression42 YM = Yi = AY 123 C---- N Number of input items -43 XM2 = XI*XT > AX*AX . 124 C---- x Independent variable44 SKX = SKX + XM*KM ‘ 125 C-+-- Y Dependent variable45 SXY = SXY + XM#YM 126 C---- A y- intercept46 SYY = SYY + YN*YM 127 C---- 8 Slope -47 SXX2 © SXXZ + XMHXMZ 128 C---- R Correlation coefficient48 SX2X2 = SK2X2 + XMZ*KMZ 129 C---- SD Standard deviation of residuals49°20. SYX2 = SYX2 + YMAXM2 130 C----A \ 131 C---- If ITYPE = 1 then linear (no prior x,y conversion}51 B= (SXY*SX2X2 - SYX2*SKX2Z) / (SXX#SK2K2 = SXX2*SXX2Z) 132 C---- = 2 then log (x=log x,y=l0g y)52 C = (SKX*SYX2 - SXX2*SKY) / (SXX*SX2X2 - SXX2*SXX2) - 133 € 7 = 3 then exponential (x=log x) -53 A= AY = BAAX = CFAX*AX s 134o4 - 135 DIMENSTON X(2), Y(2)55 C---- Standard deviation of residuals 136A aoe 137 DOUBLE PRECISION SUMX2,SUMY2, XY87 00 150 J = 1,N 13858 TY = YW) 13959 T= x(d) 14060 RV = A+ BATX + CHTK*TX - 14161 ERR = TY - RV 14262 AVGY=AYGY+TY - 143~ 63 150 S = S + ERR*ERR ~ 14464 ~ 145 00 100 I=1,N65 S = SCRT(S/(N-1)) - - 146 TX=X(1)66 . < 147 LFCLTYPE.£0.2)TX*ALOG(TX)67 RETURN - - 148 TY=¥(1)68 END 149 LF(LTYPE.NE.1)TY*ALOG(TY)= SUMK+1 STITLE: ‘REGVAL: Get regression vatue’ io SUEKDSSUNK2 STRATE2 SSUBTITLE: ‘Copyright 1986 PJ Kaufman’ 182 SUMYCSUMY4TY3 SSTORASE - 153 SUMY 2=SUMY2+TY*TY -: FUNCTION REGVAL(K,A,B,ITYPE) - 154 KY@XY4TXETY, . ~
155 100 CONTINUE
 640 APPENDIX 2 APPENDIX 2 bat
B40 RCS eet
 
 
 
186 . X LINEAR = CURVI. LOG EXPON157 LFCN.LT.2)RETURN 1.00 “61 63.54 -84158 B= (N*XY - SUMX*SUMY ) / (N*SUMX2- SUMX*SUMX ) 1.50 178 280.75 192159 Am ( SUHY -B*SUMX) /N 2.00 295, 96.98 1.02160 2.50 1.12 3.13 1.12 1.12161 C---- Correlation coefficient 3.00 1.29 1.29 1.30 1.24162 R1=N*XY = SUMX*SUMY 3.50 1.46 1.46 1.47 1.36163 R= ( N*SUMX2 - SUMK*SUMX) * (N*SUMY 2 - SUMY*SUMY ) 4.00 1.63 1.63 1.64 1.50164 R= ((RI**2)/R2)*100 4.50 1.80 1.80 1.80 1.66165 5.00 1.97 1.97 1.95 1.83166 C---- Standard deviation around fitted function 5.50 2.14 2.14 2.1 2.01167 sb=0 6.00 2.31 2.31 2.26 2.22\ 168 00 150 J=1,N 6.50 2.48 2.48 2.41 2.45169 Ty=¥(0) 7.00 2.65 2.66 2.56 2.70170 TX=K(0)17a RV=@REGVAL(TX,A,8, ITYPE)172 TeTY-RV173 150 SO=SO+T*T .174 SD=SORT(SO/(N-1))175176 RETURN77 END} LEAST-SQUARES SOLUTION FOR CORN VERSUS SOYBEANSig Least-Squares Regression Analysis Data Input‘ Y (dep) X (ind): 1.27 2.431.19 2.261.10 2.151.10 2.071.05 2.031.00 2.4598 2.361.09 2.4400 | .1.12 2.521.18 2.74 ‘1.16 2.981.24 2.931.03 2.691.08 2.631.15 2.631.33 3.081.08 3.241.57 6.222.55 6.123.02 6.332.54 4.92 ~ —2.25 6.81 - -2.02 5.88 “2.25 6.61 -2.52 6.283.11 7.522.50 6.05,r2 St devLINEAR Y= 7.274 + -339*X Q4.2 1278 - -CURVE Y= 310 + .3234X + © .002*x2 +278LOG Ln(Yo= 616 + -79S*L (XD 86.3 .282
EXPON La(Y)= 2.373 + -195*K 87.3 .283
  542 APPENDIX 2 APPENDIX 2 -
643nn
ee
——
Plot of fmoression Analysis Plot of Regression Analvsis
conn3.073.032992.952942.872ies21792762072Bee216s22602:56sis2.48Bae2402382.322282262.20216zz20820s2iot1971293189iesie.17717s1169alesignLs71531149,3.45alaAls?13a1230ze4122aeeary0tos102
98
 
   
 
 
 
    
 
 
 
  
1102 1xee x19 tx
tt
2.75 384 4.34 B14 5.94 6.73 7.5~ ‘SOYBEANS
 
 
 
 
EXPON A= ~.373. B= 119. c= +000
   
644
APPENDIX 2
SE
Piet of Kearession Anal ve:
 
 
 
 
APPENDIX 2 -
645
 
LEAST-SQUARES SOLUTION FOR SOYBEANS ONLY
Least Squares Regression AnalysisY (dep) X ¢indd
2.43 1.002.26 2.002.15 3.002.07 4.002.03 5.002.45 6.002.36 7.002.44 8.002.52 9.002.74 10.002.98 11.002.93 12.002.69 13.002.63 14.002.63 15.003.08 16.003.24 17.006.22 18.006.12 19.006.33 20.004.92 21.006.81 22.005.88 23.006.61 24.006.28 25.007.61 26.006.05 27.00reLINEAR Y= 987 + 211K 78.3CURVE Y= 1,684 + = .1124X + ,004*K2LOG Ln Y= .271 + .416*Ln(X) 58.0EXPON Ln(Y)= .527 + -053*X 83.1x LINEAR CURVI. LOG EXPON1.00 1.20 1.80 1.31 4.792.00 1.41 1.92 1.75 1.883.00 1.62 2.05 2.07 1.994.00 1.83 2.19 2.34 2.095.00 2.04 2.33 2.56 2.216.00 2.25 2.48 2.77 2.337.00 2.47 2.64 - 2.95 2.458.00 2.68 2.80 3.12 2.59- 9.00 2.89 2.98 3.27 2.7310.00 3.10 3.16 3.42 2.8711.00 3.31 3.34 3.56 9 3.0312.00 3.52 3.53 3.69 3.1913.09 3.73 3.74 3.82 3.3714,00 3.94 3.94 3.93 3.55.15.00 4.15 4.16 4.05 3.7416.00 4.37 4:38 4.16 3.9417.00 4.58 4.61 4.27 4,1618.00 4.79 4.84 4.37 4,3819.00 $.00 5.09 4.82 4,6220.00 5.21 5.34 4.56 4.8721.00 5.42 5.59 4.66 5.14
Data
Input
St dev-881830
1.212806 
646
APPENDIX 2
eA
22.23.24,25.26.27.28.23.30.31.32.33.34.35.36.37.38.39.40.
0000000000000000000000000000oaco0009000
WOOD OMOVIYIN a aMADED
-6384-05274869+90cil32537495<i?3859
80
0122243
5.866.136.41
6.99
8.238.568.899.249.599.9410.3110.6812.0611.4411.83
OOMHHHTAATAAATHN RRS
7584+92OL091722533-4024885626976-B390960309
wLMOVINaAaaaE
10.75ql.95
ul
12.13.14.
4171023469+0543-83267118
6820
3460
2800
APPENDIX 2
647
 
Plot of Keoression analvels
SUYBEANS
7.81
22407.307.207.09
a9 t
aay6.786.68.esseeas
e.37 t
e.2/e.16a.08
$196 5
ses
S273 43.65 ¢5.54}
5.44
5.34 1
3.28
5.13
S.054.924182are
Alen:
alsaalan4.504.20asiaS993.89Bury3.68sis3.483.373.27sity3.062962.8627s2652332442isa
2.24 1
21s2.03
 
 
't‘'tt't''t
 
 
 
 
  
 
  
 
 648 APPENDIX 2 APPENDIX 2 - a9i ed
————— —
                     
Plot of Kecression analvers
   
: SOYBEANS,
7.31
7.40
250: 7220
“098.99e1898.78eee8.58ea72.376.27ents
1t
 
   
   
   
650 APPENDIX 2eee EEUNIK GD
——
Plot of Reoression Analysis
 
  
: xe 4: - :t x ax: x .: x .t “xx.i
 
19.20 22.91 26.6
 
0048
 
- APPENDIX 3ba
Matrix Solution to LinearEquations and Markov Chains
DIRECT SOLUTION AND CONVERGENCE METHOD
Before computer programs offered ready solutions, problems such as Markov chains weresolved in a direct manner, by algebraically manipulating the equations. This direct solutionrequires an understanding of simple matrix arithmetic, and very careful attention to calcu-lating the numbers correctly. The convergence method is now easier, although it requiresmany more calculations. Without the computer, we would never even consider using thisapproach—with a computer, it is the best choice.
GENERAL MATRIX FORM
A matrix is a rectangular arrangement of elements into rows and columns. A matrix A issaid to be m x n (pronounced “m by n”) if there are m rows and n columns in A.
 
Qn a2 ain
x An GamGmn =| OR SP t
Bmi Gm Amn
Certain properties of a matrix make it a valuable tool for solving simultaneous linear equa-tions. These elementary matrix operations, called transformations, allow you to alter therows (which will represent equations) without changing the solution. There are threebasic row operations:1. Multiplication or division of all elements of the row by any number.2. Interchanging of any two rows (and consequently of all rows).3. .The addition or subtraction of the elements of one row with the correspondingelements of another.‘To relate the matrix to simultaneous linear equations, consider a three-equation exam-ple, where only the coefficients, a, remain as unknowns:GX t+ 2+ 4y3%s= Ay (1) —Gay + axX2 + A23%3= Gry (2)GX + As2X2 + Ass%s= Ass (3) -The three elementary row operations can now be interpreted in terms of simple oper-ations on an equation:1. When both sides of an equation are multiplied or divided by the same value, the re-sults are equal.
2. Any two equations in a system of equations can be interchanged with no effect.3. When two equals are added or subtracted, the results are equal.
DIRECT SOLUTICN
Putting these rules into use, write the coefficients of the three simultaneous linear equa-tions as a 3 x 4 coefficient matrix, 
 
APPENDIX 3
@y G2 43 Gu\ ~@n Gr 13 au 13, @s2 As3 Gs
with the objective of reducing the matrix to the form
100 A01040 0 1 As
which would mean that x; = Ay, %2, = 42, and x5 = As since
652
1X2x,+0X 22+ 0% %3=A1Ox 2x + 1X 22+ 0 X23 = Ad
Ox x, + OX a+ 1553
‘To achieve the results, divide the first equation by 41, the first element, leaving
(: an as a) an
4, 4n an
and then multiply by an to get the first elements in rows one and two the same:
Q2X@n @3XGn aXG(on 1 2) aayan au au
Now subtract (1-2) from (2) and get
anXan ayX an au Xan(0 Gq — EA ay a (2-1)an - au au
‘That successfully eliminates the first element az; from the second equation. By going backand multiplying (1,1) by as: and subtracting the resulting equation from (3), 43; can beeliminated from equation (3). Now column 1 looks like
() Co
Column two or any column can be operated upon in the same manner as the first column:
1. Divide row n by the element in position n (for-example, @nn), thereby settingAnn = 1.
2. Maltiply row n by the corresponding element in row 7 # 7, SO that don X Gin = Fine
3. Subtract row 1 from row 3, resulting in a, = 0, and all other elements reduced by
the corresponding element in row 7.Continue this procedure