10 CLS:PRINT "MULTFIT-Fits Multiplicative Model"
11 PRINT" by Michael Birnbaum"
12 PRINT :PRINT
14 PRINT"Prediction(i,j) = B + S(i)T(j)"
15 PRINT :PRINT "Initial Estimates as Follows"
16 PRINT" B estimated from 4 corners"
17 PRINT " S and T estimated from marginal sums"
18 PRINT
20 INPUT"NO. OF ROWS";r
30 INPUT"NO. OF COLS";c
40 DIM x(r,c),s(r),T(c),S1(r),T1(c),S2(r),T2(c),X2(r,c)
41 INPUT"Data in file (f) or Entered from keyboard (k)";k$
42 IF k$ <> "f" AND k$ <> "k" THEN PRINT "type either 'k' or 'f': go to 42
43 IF k$ = "f" OR k$ = "F" THEN 1000
45 PRINT :PRINT "ENTER DATA, BY ROWS":PRINT
50 FOR i=1 TO r
55 PRINT :PRINT "ROW";i
60 FOR j=1 TO c
70 PRINT "X(";i;",";j;")=";
80 INPUT x(i,j)
90 NEXT:NEXT
120 FOR i=1 TO r:S1(i)=0:S2(i)=0 :NEXT
125 FOR j=1 TO c:T1(j)=0:T2(j)=0:NEXT
130 A1=(x(r,c)-x(r,1))/(x(1,c)-x(1,1))
132 IF A1=1 THEN A1=1.0000001#
135 B=(x(1,1)*A1-x(r,1))/(A1-1)
140 CLS:PRINT "FIRST EST. OF B = "
141 PRINT B:PRINT
150 FOR i=1 TO r:FOR j=1 TO c:X2(i,j)=x(i,j)-B:NEXT:NEXT
155 GT=0
160 FOR i=1 TO r:S1(i)=0:FOR j=1 TO c
170 GT=GT+X2(i,j)
175 S1(i)=S1(i)+X2(i,j)
180 T1(j)=T1(j)+X2(i,j)
185 NEXT:NEXT
190 IF GT=0 THEN PRINT"Grand Total=0. Reflect Data":END
191 PRINT :PRINT "ESTIMATES OF S"
192 FOR i=1 TO r:s(i)=S1(i)/GT:PRINT i,s(i):NEXT
193 PRINT "HIT A KEY"
194 a$=INKEY$:IF a$="" THEN 194
195 PRINT :PRINT "ESTIMATES OF T":PRINT
196 FOR j=1 TO c:T(j)=T1(j):PRINT j,INT(T(j)*1000+.5)/1000:NEXT
198 PRINT: PRINT "EST. OF B =";B: PRINT
203 PRINT "HIT A KEY"
204 a$=INKEY$:IF a$=""THEN 204
300 ICT=0:SER=0
310 CLS:PRINT "FIT OF INITIAL ESTS"
315 PRINT " I J ","X(I,J)"," PRED"," RES"
320 FOR i=1 TO r:FOR j=1 TO c: ICT = ICT+1
322 PR=s(i)*T(j)+B
324 E1=x(i,j)-PR:SER=SER+E1*E1
326 PRINT i;j,x(i,j),INT(1000*PR+.5)/1000,INT(1000*E1+.5)/1000
327 IF 14 > ICT THEN 330
328 ICT = 0: PRINT "HIT A KEY"
329 a$=INKEY$:IF a$=""THEN 329
330 NEXT:NEXT
332 PRINT "HIT A KEY";
335 a$=INKEY$:IF a$=""THEN 335
340 PRINT :PRINT "SUM OF SQUARED ERRORS":PRINT "SS=";SER
347 INPUT"ANOTHER LOOK (Y/N)?";r$
348 IF r$="Y" OR r$="y" THEN GOTO 191
350 CLS: PRINT"Iterative Least-Squares Solution"
352 INPUT"Tolerance (e.g., .01)"; TL
360 REM -----another iteration starts here ------
365 BT = 0
370 FOR i = 1 TO r: FOR j=1 TO c: BT=BT+x(i,j)-s(i)*T(j):NEXT:NEXT
375 B=BT/(r*c)
380 FOR i = 1 TO r: S1(i)=0:S2(i)=0: FOR j=1 TO c
382 S1(i)=S1(i)+(x(i,j)-B)*T(j)
383 S2(i)=S2(i)+T(j)*T(j)
384 NEXT: S1(i)=S1(i)/S2(i):NEXT
390 FOR j=1 TO c: T1(j)=0:T2(j)=0:FOR i=1 TO r
392 T1(j) = T1(j)+(x(i,j)-B)*s(i)
393 T2(j) = T2(j)+s(i)*s(i)
394 NEXT: T1(j)=T1(j)/T2(j):NEXT
400 SD = 0: SER =0
410 FOR i=1 TO r: SD=SD+ABS(S1(i)-s(i)):NEXT
420 FOR j=1 TO c: SD=SD+ABS(T1(j)-T(j)):NEXT
430 FOR i=1 TO r: FOR j=1 TO c
432 SER=SER+(x(i,j)-B-s(i)*T(j))^2
435 NEXT:NEXT
440 PRINT"Change in Scales=";SD, "Sum of Squares=";SER
450 FOR i=1 TO r: s(i)=S1(i):NEXT
460 FOR j=1 TO c:T(j)=T1(j):NEXT
465 IF TL > SD THEN GOSUB 700
470 IF SD > TL THEN 360
482 INPUT"Another look (y/n)";RL$
484 IF RL$="Y" OR RL$="y" THEN GOSUB 700
486 INPUT"More iterations (y/n)";AM$
488 IF AM$="Y" OR AM$="y" THEN 350
490 INPUT"Save results (y/n)";AD$
498 IF AD$="Y" OR AD$="y" THEN GOSUB 900
500 END
700 PRINT :INPUT"HIT A KEY";KK$
710 CLS: PRINT"Approximate Least-Squares Solution"
791 PRINT :PRINT "ESTIMATES OF S"
792 FOR i=1 TO r: PRINT i,INT(s(i)*1000+.5)/1000:NEXT
793 PRINT "HIT A KEY"
794 a$=INKEY$:IF a$="" THEN 794
795 PRINT :PRINT "ESTIMATES OF T":PRINT
796 FOR j=1 TO c: PRINT j,INT(T(j)*1000+.5)/1000:NEXT
798 PRINT: PRINT "EST. OF B =";B: PRINT
803 PRINT "HIT A KEY"
804 a$=INKEY$:IF a$=""THEN 804
810 CLS:PRINT "Approximate Least Squares Fit":PRINT
812 SER=0
815 PRINT " I J ","X(I,J)"," PRED"," RES"
820 FOR i=1 TO r:FOR j=1 TO c
824 PR=s(i)*T(j)+B
826 E1=x(i,j)-PR:SER=SER+E1*E1
828 PRINT i;j,x(i,j),INT(1000*PR+.5)/1000,INT(1000*E1+.5)/1000
830 NEXT
832 PRINT "HIT A KEY"
835 a$=INKEY$:IF a$=""THEN 835
838 NEXT
840 PRINT :PRINT "SUM OF SQUARED ERRORS":PRINT "SS=";SER
850 RETURN
900 REM subroutine to create file
901 INPUT"type 'y' if you want to save output (y/n)";pout$
902 IF pout$ <> "y" AND pout$ <>"Y" THEN GOTO 999
903 INPUT"name of file";nemfile$
904 INPUT"title of output";title$
908 OPEN nemfile$ FOR OUTPUT AS #1
910 PRINT #1, "Approximate Least-Squares Solution Multfit"
912 PRINT #1, title$
914 PRINT#1, :PRINT #1, "ESTIMATES OF S"
916 FOR i=1 TO r: PRINT #1, i,INT(s(i)*1000+.5)/1000:NEXT
920 PRINT #1,:PRINT #1,"ESTIMATES OF T":PRINT #1,
922 FOR j=1 TO c: PRINT#1, j,INT(T(j)*1000+.5)/1000:NEXT
930 PRINT# 1,: PRINT #1, "EST. OF B =";B: PRINT #1,
932 PRINT # 1, "Approximate Least Squares Fit":PRINT#1,
934 SER=0
936 PRINT #1, " I J ","X(I,J)"," PRED"," RES"
938 FOR i=1 TO r:FOR j=1 TO c
940 PR=s(i)*T(j)+B
942 E1=x(i,j)-PR:SER=SER+E1*E1
944 PRINT #1,i;j,x(i,j),INT(1000*PR+.5)/1000,INT(1000*E1+.5)/1000
946 NEXT
948 NEXT
950 PRINT#1, :PRINT#1, "SUM OF SQUARED ERRORS":PRINT #1, "SS=";SER
955 PRINT#1,
960 PRINT#1, "Extra Copy of Predictions (Copy and Paste)"
970 FOR i=1 TO r:FOR j=1 TO c
974 PR=s(i)*T(j)+B
978 PRINT #1,INT(1000*PR+.5)/1000
980 NEXT
982 NEXT
983 PRINT#1,
985 PRINT#1, "Extra Copy of DATA (Copy and Paste)"
987 FOR i=1 TO r:FOR j=1 TO c
990 PRINT #1,x(i,j)
992 NEXT
993 NEXT
995 CLOSE #1
997 RETURN
999 END
1000 INPUT"Your data are in a file (y/n)";a$
1001 IF a$ <> "y" AND a$ <> "Y" THEN 45
1020 INPUT"Name of file ="; dmul$
1023 OPEN dmul$ FOR INPUT AS #2
1026 FOR i = 1 TO r
1028 FOR j = 1 TO c
1030 INPUT #2, x(i,j)
1034 NEXT:NEXT
1040 CLOSE #2
1050 GOTO 120