1  ! PROGRAM FQEQ, VERSION 1A (31-AUG-73)
   !
2  ! PROGRAM FQEQ IS A PROGRAM IN THE PROBIT ANALYSIS CHAIN.
   ! IT IS ONE OF THE SEVERAL PROGRAMS AND FILES THAT COMPOSE COSAP,
   ! LAWRENCE UNIVERSITY'S CONVERSATIONALLY ORIENTED STATISTICAL
3  ! ANALYSIS PACKAGE.  COSAP IS DISTRIBUTED FOR EDUCATIONAL USE, ON
   ! THE CONDITION THAT IT NOT BE SOLD, RENTED OR LEASED FOR PROFIT.
   !
4  ! IF YOU DESIRE FURTHER INFORMATION ABOUT COSAP, PLEASE CONTACT THE
   ! LAWRENCE UNIVERSITY COMPUTER CENTER, APPLETON, WISCONSIN, 54911.
   !
10 DIM N%(30)  ! FUNCTION TO CREATE FILE EXT. BASED ON JOB#
25 DEF FNZ$ : IF N%(0)=30% THEN 40 
    ELSE CHANGE SYS(CHR$(6)+CHR$(9)+CHR$(0)) TO N% : 
   Z8$=".J"+RIGHT(NUM$(N%(1)/2%),2%)
40 FNZ$=Z8$! EXTENSION IS".J"+JOB #
55 FNEND
70 !****FQEQ
85 OPEN "DAT1"+FNZ$ AS FILE 1
100 DIM#1,X(200,50),H$(30)=2,C(50,50),M4(0),M5(0),M9(0),N(0),J9(0),
          I4(0),I7(0),N6(0),X2(50),L1$(8)=32%
115 OPEN "PROB1"+FNZ$ AS FILE 2
130 DIM#2,Y2(200),W1(200),W2(200),A(0),B(0),V3(0),S7(0),T7(0),R7(0)
145 S7=S7(0):T7=T7(0):R7=R7(0)
160 DEF FNR(Q)
175 INPUT Z$
190 Q=INSTR(1,Z$,"SOS")
205 IF Q>0 THEN 250 
220 Q=VAL(Z$)
235 GOTO 265 
250 Q=808
265 FNR=Q
280 FNEND
295 V3=0:PRINT
310 PRINT"DO YOU WISH TO SPECIFY AN ESTIMATE OF THE NATURAL RESPONCE RATE C";
325 INPUT N$:IF N$="YES" THEN 400 
340 IFN$<>"SOS" THEN 460 
355 PRINT:PRINT"C IS THE NATURAL MORTALITY RATE AND WILL BE USED TO"
370 PRINT"CONTROL FOR ANY RESPONSES WHICH MIGHT HAVE OCCURED"
385 PRINT"WITHOUT THE STIMULI.":GOTO 310 
400 PRINT:PRINT"SPECIFY THE VALUE OF C (MUST BE LESS THAN 1 !)";
415 V3=FNR(Q):IF V3<>N6(0) THEN 445 
430 PRINT:PRINT"IF YOU ARE UNCERTAIN, SPECIFY 0.":GOTO 400
445 IF 1<=V3 THEN 400 
460 W9,X3,S3,S2,S1=0
475 FOR I=1 TO N(0):P=W1(I)
490 IF P<=0 THEN 685 
505 IF P>=1 THEN 685 
520 W9=W9+1
535 D=P
550 IF D<=.5 THEN 580 
565 D=1-D
580 T2=LOG(1/(D*D))
595 T=SQR(T2)
610 A2=2.515517+.802853*T+.010328*T2
625 D2=1+1.432788*T+.189269*T2+.001308*T*T2
640 Z=T-A2/D2
655 IF P>.5 THEN 670 ELSE Z=-Z
670 Z=Z+5:X3=X3+Y2(I):S3=S3+Z:LET S2=S2+Y2(I)^2:S1=S1+Y2(I)*Z
685 NEXT I
700 B=(S1-(X3*S3)/W9)/(S2-(X3*X3)/W9)
715 X3=X3/W9:S3=S3/W9:A=S3-B*X3:LET D2=0
730 FOR I=1 TO N(0):W2(I)=A+B*Y2(I):NEXT I
745 W8,W5,W3,W4,W6=0
760 FOR I=1 TO N(0):Y=W2(I):D=Y-5:LET A9=ABS(D)
775 T=1/(1+.2316419*A9)
790 Z=.3989423*EXP(-D*D/2)
805 T1=.3193815*T:T2=.3565638*T^2:T3=1.781478*T^3
820 T4=1.821256*T^4:T5=1.330274*T^5
835 P=1-Z*(T1-T2+T3-T4+T5)
850 IF D>=0 THEN 880 
865 P=1-P
880 Q1=1-P
895 W=(Z*Z)/(Q1*(P+V3/(1-V3)))
910 IF Y>5 THEN 955 
925 P3=(Y-P/Z)+W1(I)/Z
940 GOTO 970 
955 P3=(Y+Q1/Z)-(1-W1(I))/Z
970 W9=W*X(I,T7)
985 W8=W8+W9:W5=W5+W9*Y2(I):W3=W3+W9*P3
1000 W4=W4+W9*Y2(I)^2:W6=W6+W9*Y2(I)*P3
1015 NEXT I
1030 X3=W5/W8:S2=W4-W5*W5/W8:S1=W6-W5*W3/W8
1045 B=S1/S2:A=W3/W8-B*X3:S2=0
1060 FOR I=1 TO N(0)
1075 Y=A+B*Y2(I)
1090 D=W2(I)-Y
1105 S2=S2+D*D
1120 W2(I)=Y
1135 NEXT I
1150 IF ABS(D2-S2)<=1E-7 THEN 1180 
1165 D2=S2:GOTO 745 
1180 A(0)=A:B(0)=B:V3(0)=V3
1195 CLOSE 1,2
1210 K$=SYS(CHR$(7)):N$=LEFT(K$,8):S0=VAL(RIGHT(K$,9))
1225 CHAIN N$ S0
1240 END