Trailing-Edge
-
PDP-10 Archives
-
decuslib20-03
-
decus/20-0093/shgm.for
There are no other files named shgm.for in the archive.
C MODEL SUBROUTINE FOR SCHIZ SIMULATION
C OCTOBER 19, 1971
C BOB STOUT, PROGRAMMER
C D. MALIN, SIMULATION DESIGNER
C MODIFIED AUGUST, 1972 BY BOB STOUT
C MODIFIED JANUARY, 1973 BY BOB STOUT
C
C
SUBROUTINE MODEL(N,NCOND,*)
C
INTEGER*4 SI(12),SF(12), OVNAM
INTEGER*4 IV(12), SINK
LOGICAL*4 COSTPT, DATOUT
REAL*4 FV(12)
REAL*4 EXCOST(9),ECOSTF(2),SCOSTB(9),SCOSTI(2),MCOSTF(2),ACOSTF(3)
REAL*4 OHDRT(3),SREIM,TCOST(2),TOCST(2)
INTEGER*4 YES1,YESY,ERSP
!! DATA YES1/'1'/, YESY/'Y'/
INTEGER*4 INST,SAMP,NT,EXTER,PCOST,IX,IR
LOGICAL*4 LEAVE,SKIZ,MAXPT,PNTIT
COMMON /VARS1/IV,FV,SI,SF
COMMON /IO1/NOV,SINK,OVNAM(8,6),COSTPT,DATOUT
COMMON /IO/IDEV1,IDEV2,IDEV3,IDEV4
C
INTEGER*4 REL,ADOPT,DIAG,SEXI,SEXS
EQUIVALENCE (INDEX,IV(1)),(REL,IV(2)),(ADOPT,IV(3)),
1 (DIAG,IV(4)),(SEXI,IV(5)),(SEXS,IV(6)),(METHOD,IV(7)),
2 (INST,IV(8)),(SAMP,IV(9)),(EXTER,IV(10)),(PCOST,IV(11))
C
REAL*4 METH(2), NOISE
C METH CONTAINS THE MAIN EFFECT PARAMETERS FOR METHOD
C BASE IS THE BASE RATE OF SCHIZ IN THE GENERAL POPULATION
C DMINT IS THE COEFFICIENT FOR THE DIAG-METHOD INTERACTION
C
C SEXEFF CONTAINS PARAMETERS FOR THE SEX EFFECT. EFFECT DEPENDS
C ONLY ON WHETHER SEXI(INDEX S)=SEXS(ACTUAL S). FOR SEXI=SEXS
C SEXEFF(2) IS USED, OTHERWISE SEXEFF(1).
REAL*4 SEXEFF(2)
C
C RATES CONTAINS THE SCHIZ RATES FOR VARIOUS RELATIVES: 1=UNREL,
C 2=AUNTUNCL, 3=COUSIN, 4=SIBLING, 5=FRATTWIN, 6=IDTWIN, 7=PARENT,
C 8=CHILD, 9=NIECNEPH
REAL*4 RATES(9)
C
C DIAGEF AND ADOEFF CONTAIN THE MAIN EFFECT PARAMETERS FOR STRICY-
C NESS OF DIAGNOSIS AND ADOPT
REAL*4 DIAGEF(2), ADOEFF(3)
C
C------------------------------
C
DOUBLE PRECISION MODNM1,MODNM2,MODNM3
COMMON /DATANM/ MODNM1,MODNM2,MODNM3
DATA MODNM1/'SHGM.DAT'/,MODNM2/'SHGM3.DAT'/,MODNM3/'SHGM.BIN'/
DATA YES1/'1'/, YESY/'Y'/
METHOD=IV(7)+SI(7)
TOHRT=OHDRT(INST)+TOCST(SAMP)
NOISE=NOISE!DEC-10 REQUIES THIS
NCOND=NCOND!DEC-10 BUG MAKES THIS NECESSARY
XN=N
SEX=0.
C COMPUTE OVERALL SEX EFFECT
IF(SEXI.EQ.0 .OR. SEXS.EQ.0)GO TO 2
C
C IT ONLY MATTERS WHETHER SEXI=SEXS, NOT WHAT THE SEXES ARE
SEX=SEXEFF(1)
IF(SEXI.EQ.SEXS)SEX=SEXEFF(2)
C IDENTICAL TWINS ALWAYS HAVE SAME SEX
2 IF(REL.EQ.6)SEX=SEXEFF(2)
IF(SEX.GT.0.01)GO TO 3
C
C SEX IS RANDOMIZED--COMPUTE AVERAGE EFFECT
CALL BINOM(N,.5,I)
XI=I
SEX=XI*SEXEFF(1)+(XN-XI)*SEXEFF(2)
SEX=SEX/XN
C
C------------------------------
C
C ***** THE MODEL *****
C
C FIRST SET UP SOME THINGS
C
3 LEAVE=.FALSE.
P1=BASE*(1.+DIAGEF(DIAG))
C
C COMPUTE BASIC COST FACTORS
C
ECST=EXCOST(REL)+ECOSTF(INDEX)+TCOST(SAMP)
SCST1=SCOSTB(REL)+SCOSTI(INDEX)
SCST2=SCST1+MCOSTF(METHOD)*SCST1
SCST3=SCST1+ACOSTF(ADOPT)*SCST1
SCOST=SCST1+SCST2+SCST3+SREIM
C SKIP THE MODEL IF INDEX=NONSCHIZ OR REL=UNREL
IF(INDEX.EQ.2 .OR. REL.EQ.1)GO TO 4
C
C HEREDITY-REARING FACTORS
P1=RATES(REL)*ADOEFF(ADOPT)
C
C SEX LINKAGE FACTOR
P1=P1+(SEX-1.)*P1*(1.-P1)
C
C ADD IN BASE RATE
P1=P1+BASE*(1.-P1)
C
C PUT IN GROSS RATE INFLATION FACTORS--METHOD AND DIAGNOSIS
P2=1.-P1
P2=P1*P2*P2
P1=P1+(METH(METHOD)+DIAGEF(DIAG)+DMINT*DIAGEF(DIAG)*METH(METHOD))
1 *P2
C
C SELECT SAMPLING FUNCTION
C
4 GO TO(600,500,700),SAMP
C
C NO. OF SUBJECTS DIAGNOSED SCHIZ IS BINOMIAL RV
600 LEAVE=.TRUE.
NT=N
CALL BINOM(N,P1,NSCHIZ)
C
C PRINT RESULTS
WRITE(SINK,1001)NSCHIZ
GO TO 550
C
C PASCAL SAMPLING IS CHOSEN
700 LEAVE=.TRUE.
C
C IR= # DESIRED SCHIZ; Q=PR(NONSCHIZ)
C
IR=N
Q=1.-P1
CALL PASCAL(IR,Q,IX,NT)
C
C PRINT IT ALL
C
WRITE(SINK,3020) IR,NT
3020 FORMAT(/' ',I8,' SCHIZOPHRENICS WERE FOUND UNDER EXPERIMENTER.'
1 ' ',I8,' SUBJECT-PAIRS (N) WERE SAMPLED TO DO SO.')
C
C PRINT COSTS
C
550 TSCST=SCOST*NT
OHCST=(ECST+TSCST)*TOHRT
TCST=ECST+TSCST+OHCST
WRITE(SINK,3001)
3001 FORMAT(//' COSTS:'/)
PNTIT=LEAVE .OR. MAXPT
IF(PNTIT) WRITE(SINK,3002) ECST
3002 FORMAT(' EXPERIMENTAL CONDITION SET-UP COST:',F10.2)
IF(PNTIT) WRITE(SINK,3003) NT,TSCST
3003 FORMAT(' SAMPLED N= ',I5,' YIELDS SUBJECT COSTS OF: ',F10.2)
IF(PNTIT) WRITE(SINK,3004) OHCST
3004 FORMAT(' OVERHEAD COSTS ARE: ',F10.2)
WRITE(SINK,3005) TCST
3005 FORMAT(' THE TOTAL COST IS: ',F10.2)
IF(LEAVE) RETURN
GO TO 502
C
1001 FORMAT(1H0,5X,I4,' SUBJECTS DIAGNOSED SCHIZOPHRENIC.')
C
C
C NOW FOR THE BERNOULLI-ONE-AT-A-TIME SAMPLING
C
500 NT=0
NFS=0
P9=1.-P1
IF(P9.LT.0) P9=0
IF(P9.GT.1) P9=1
MAXPT=EXTER.EQ.1 .OR. PCOST .EQ.1
IF(EXTER.EQ.1) WRITE(SINK,3011)
3011 FORMAT(//' WELCOME TO THE SCHIZGAME SEQUENCE WHICH ALLOWS'/
1' YOU, THE EXPERIMENTER, TO DECIDE TO CONTINUE'/
2' EXPERIMENTATION OR STOP. AT THE END OF EACH'/
3' SUBJECT-PAIR SAMPLED, YOU WILL BE ASKED IF YOU'/
4' WISH TO CONTINUE THE EXPERIMENT OR NOT. IF YOU'/
5' DECIDE TO CONTINUE, ENTER A "1"; IF NOT, A "0".'//)
IF(PCOST.NE.3 .OR. EXTER.EQ.1) WRITE(SINK,3017) ECST,SCOST,TOHRT
3017 FORMAT(//' BASE COSTS:'/' EXPERIMENTAL CONDITION SET-UP: ',F10.2/
1 ' PER SUBJECT-PAIR SAMPLED COST: ',F10.2/' OVERHEAD RATE: ',F5.3)
GO TO 501
502 WRITE(SINK,3012)
3012 FORMAT(//' CONTINUE?--"1"=YES, OR "0"=NO.')
READ(IDEV1,3013)ERSP
3013 FORMAT(A1)
LEAVE=ERSP.NE.YES1 .AND. ERSP.NE.YESY
IF(LEAVE) GO TO 550
501 NT=NT+1
SKIZ=.FALSE.
CALL CHANCE(P9,&505)
SKIZ=.TRUE.
505 IF(SKIZ) GO TO 506
WRITE(SINK,3014) NT
3014 FORMAT(/' SAMPLED SUBJECT-PAIR NUMBER: ',I5/
1 ' THE ACTUAL SUBJECT WAS NOT FOUND TO BE SCHIZOPHRENIC.')
507 WRITE(SINK,3015) NFS,NT
3015 FORMAT(/' NUMBER OF SCHIZOPHRENICS FOUND SO FAR: ',I5/
1 ' TOTAL NUMBER OF SUBJECT-PAIRS SAMPLED SO FAR: ',I5)
GO TO 550
506 WRITE(SINK,3016) NT
NFS=NFS+1
3016 FORMAT(/' SAMPLED SUBJECT-PAIR NUMBER: ',I5/
1 ' THE ACTUAL SUBJECT IS A SCHIZOPHRENIC.')
WRITE(SINK,3015) NFS,NT
GO TO 550
C
ENTRY MINIT
OPEN(UNIT=IDEV3,ACCESS='SEQIN',FILE=MODNM2,DIRECTORY='4030,11')
C
READ(IDEV3,2000) EXCOST,ECOSTF,SCOSTB,SCOSTI,MCOSTF,ACOSTF,OHDRT,
1 SREIM,TCOST,TOCST,
1 METH,BASE,DMINT,SEXEFF,RATES,DIAGEF,ADOEFF
2000 FORMAT(9F8.2/2F8.2/9F8.2/2(2F8.2/),3F8.2/3F8.2/F8.2/2F8.2/2F8.2/,
1 2F8.2/2(F8.2/),2F8.2/9F8.2/2F8.2/3F8.2)
RETURN
END
C SUBROUTINE FOR THE NEGATIVE BINOMIAL OR PASCAL DISTRIBUTION
C FOR SAMPLING THE N BERNOULLI TRIALS NEED TO OBTAIN
C R (IR) SUCCESSES SPECIFIED BY THE USER. THE SUBROUTINE
C RETURNS BOTH THE VALUE FOR N (IN IN), AND THE VALUE FOR
C X (IN IX), WHICH IS THE NUMBER OF FAILURES (ZEROS) IN THE
C N TRIALS. P IS THE PR OF A FAILURE (ZERO).
C
SUBROUTINE PASCAL(IR,P,IX,IN)
INTEGER*4 IR,IX,IN
REAL*4 P,Q,B,GVAL
C
C ERROR CHECKS AND FIX-UP
C
IF(IR.LE.0) STOP 4001
IF(P .LT. -.00001 .OR. P .GT. 1.00001) STOP 4002
IF(P.LT. (.0000001)) P=.0000001
IF(P.GT. (.9999999)) P=.9999999
C
C START COMPUTATION
C
Q=1.-P
C
C MAKE SURE WE CAN NOT EXPECT OVERFLOW
C
IF(((IR*Q/P)+IR) .GT. 2000000000) STOP 4003
B=P/Q
C
C X (OR IX) IS PASCAL MEANS WE CAN USE CONVOLUTION
C OF THE GAMMA AND POISSON (MIDAS SIMULATION MANUAL)
C
CALL RGAMMA(IR,B,GVAL)
CALL POISSN(GVAL,IX)
IN=IX+IR
RETURN
END