SUBROUTINE CALCPR(NPAR, X, IERR, ICH, IALT, II, ICDAT, IR, RCDAT, 1 PROB, IUSER, RUSER, MNPCDF) C C *** THIS SUBROUTINE CALCULATES A PROBABILITY FOR THE MODEL AND *** C *** DATA GIVEN. FOR MULTINOMIAL PROBIT SOME ADDITIONAL STORAGE *** C *** CUSTOMIZATION IS REQUIRED. THIS APPROACH CAN BE *** C *** USED FOR OTHER CHOICE MODELS, WITH APPROPRIATE MODIFICATIONS *** C *** TO THE ARRAYS USED BELOW. *** C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ C INTEGER NPAR, IERR, ICH, IALT, II, ICDAT(*), IR, IUSER(*) DOUBLE PRECISION X(NPAR), RCDAT(*), PROB, RUSER(*) EXTERNAL MNPCDF C C *** CALCPR PARAMETER USAGE *** C C IALT.... NUMBER OF CHOICES AVAILABLE IN THE CHOICE SET. C ICDAT... VECTOR OF INTEGER DATA VALUES. C ICH..... INTEGER INDICATING THE CHOICE. 1 <= ICH <= IALT. C IERR.... INTEGER FOR PASSING ERROR INFORMATION. C IN THIS ROUTINE, IF IERR = 1 ON RETURN THEN THERE WERE C NO PROBLEMS. C IF IERR = 0 ON RETURN, THEN THE PROBABILITY COULD NOT C BE COMPUTED USING THE CURRENT PARAMETERS IN X. C II...... NUMBER OF INTEGER VALUES STORED IN VECTIR ICDAT. C IUSER... MODEL-RELATED INTEGER VALUES USED BY CALCPR. CONTAINS C ARRAY POINTERS TO MANAGE DATA STORAGE, AND OTHER C PARAMETERS. C MNPCDF.. SUBROUTINE WHICH CALCULATES THE CDF OF A MULTIVARIATE C NORMAL DISTRIBUTION. C NPAR.... NUMBER OF PARAMETERS IN VECTOR X. C PROB.... ON RETURN, CHOICE PROBABILITY COMPUTED USING PARAMETERS IN C X AND DATA IN ICDAT AND RCDAT. C RCDAT... VECTOR OF REAL DATA VALUES. C RUSER... MODEL-RELATED REAL VALUES USED BY CALCPR. CAN CONTAIN C USEFUL PARAMETERS, AND ALSO EXTRA WORK STORAGE. C EXTERNAL CALCP1 INTEGER ISIGU, IW, NALT, NW C ISIGU = IUSER(16) IW = IUSER(15) NALT = MAX(1,IUSER(6)) NW = MAX(1, IUSER(7)) CALL CALCP1(NPAR, X, IERR, ICH, IALT, II, ICDAT, IR, RCDAT, 1 PROB, IUSER, RUSER, NALT, RUSER(ISIGU), 2 NW, RUSER(IW), MNPCDF) C *** LAST LINE OF CALCPR FOLLOWS *** END SUBROUTINE CALCP1(NPAR, X, IERR, ICH, IALT, II, ICDAT, IR, RCDAT, 1 PROB, IUSER, RUSER, NALT, SIGU, NW, W, MNPCDF) C C *** THIS SUBROUTINE CALCULATES A PROBABILITY FOR THE MNP MODEL *** C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ C INTEGER NPAR, IERR, ICH, IALT, II, ICDAT(*), IR, IUSER(*), NALT, 1 NATTR, NW DOUBLE PRECISION X(NPAR), RCDAT(*), PROB, RUSER(*), 1 SIGU(NALT,NALT), W(NW,NALT) EXTERNAL MNPCDF C C *** CALCP1 PARAMETER USAGE *** C C IALT.... NUMBER OF CHOICES AVAILABLE IN THE CHOICE SET. C ICDAT... VECTOR OF INTEGER DATA VALUES. C IN THIS SUBROUTINE, ICDAT STORES INTEGER INDEXES WHICH C DEFINE WHICH OF THE NOMINAL ALTERNATIVES ARE AVAILABLE C IN THE CHOICE SET. (THIS IS FOR THE CASE WHEN THERE C ARE NALT NOMINAL CHOICE ALTERNATIVES, BUT NOT ALL OF C THEM NECESSARILY APPEAR IN EVERY SUBSET. IF ALL NALT C ALTERNATIVES APPEAR IN ALL SUBSETS, THEN ICSET = NALT >0 C SHOULD BE USED WITH IDUM = 1. C ICH..... INTEGER INDICATING THE CHOICE. 1 <= ICH <= IALT. C ICOV.... INDICATOR FOR TYPE OF ALTERNATIVE-SPECIFIC ERRORS, C = 0 FOR IID ERRORS, = 1 FOR CORRELATED ERRORS. C IF ICSET .NE. 0, THEN THE SAME CORRELATION MATRIX IS C USED FOR EVERY SUBSET. OTHERWISE, INTEGER DATA SHOULD C BE USED TO IDENTIFY THE ALTERNATIVES IN EACH CHOICE SET. C (STORED IN IUSER.) C IDUM... INDICATOR FOR ALTERNATIVE-SPECIFIC DUMMIES, C = 0 FOR NO, = 1 FOR YES. IF ICSET .NE. 0, THEN C THE SAME SET OF DUMMIES IS USED FOR EACH CHOICE SET. C OTHERWISE, INTEGER DATA SHOULD BE USED TO IDENTIFY THE C ALTERNATIVES IN EACH CHOICE SET (SEE NALT BELOW). C (STORED IN IUSER). C IERR.... INTEGER FOR PASSING ERROR INFORMATION. C IN THIS ROUTINE, IF IERR = 1 ON RETURN THEN THERE WERE C NO PROBLEMS. C IF IERR = 0 ON RETURN, THEN THE PROBABILITY COULD NOT C BE COMPUTED USING THE CURRENT PARAMETERS IN X. C II...... NUMBER OF INTEGER VALUES STORED IN VECTIR ICDAT. C IUSER... MODEL-RELATED INTEGER VALUES USED BY CALCPR. THE FIRST C PORTION OF IUSER CONTAINS SUCH THINGS AS ARRAY POINTERS. C IUSER ALSO CONTAINS STORED VALUES OF NATTR, IDUM, ETC. C IR...... NUMBER OF REAL VALUES STORED IN VECTOR IRDAT. C ITASTE.. INDICATOR FOR TASTE VARIATION, C = 0 FOR NO TASTE VARIATION, = 1 FOR UNCORRELATED TASTE C VARIATION, = 2 FOR CORRELATED TASTE VARIATION. C (STORED IN IUSER.) C NPAR.... NUMBER OF PARAMETERS IN VECTOR X. C PROB.... ON RETURN, CHOICE PROBABILITY COMPUTED USING PARAMETERS IN C X AND DATA IN ICDAT AND RCDAT. C RCDAT... VECTOR OF REAL DATA VALUES. C IN THIS SUBROUTINE, THE NUMBER OF DATA VALUES SHOULD C BE = IALT * NATTR SO THAT THE "GENERIC" PART OF THE C SCALE VALUE V MAY BE COMPUTED. C NALT.... TOTAL NUMBER OF NOMINAL CHOICE ALTERNATIVES (IF APPLICABLE). C IF ICSET .NE. 0 AND IDUM = 1 OR ICOV = 1 (OR BOTH), THEN C NALT SHOULD BE EQUAL TO ICSET. C OTHERWISE, NALT SHOULD BE > 0 IF EITHER IDUM OR ICOV C (OR BOTH) ARE > 0, AND ICDAT SHOULD BE USED TO PASS C INDEX INFORMATION (SEE ICDAT ABOVE). C NATTR... NUMBER OF ATTRIBUTES (I.E., REAL DATA VARS.) PER ALTERNATIVE. C NW...... NUMBER OF ROWS IN THE WORK-ARRAY W. C RUSER... MODEL-RELATED REAL VALUES USED BY CALCPR. FOR THIS MODEL, C IT CONTAINS A CONSTANT FOR THE COVARIANCE MATRIX SCALE, C AND INFORMATION USED FOR COMPUTING STEP SIZES IN FINITE- C DIFFERENCE CALCULATIONS. C SIGU.... MATRIX CONTAINING THE "UNPACKED" THE FULL COVARIANCE MATRIX C FOR ALL NALT ALTERNATIVE-SPECIFIC ERROR TERMS. THE C MATRIX IS OF DIMENSION 2 TO FACILITATE CODING. THE C NORMALIZATION USED LEAVES A ROW OF ZEROS IN THE LAST C (NALT) ROW. IT IS COMPUTED BEFORE THE CALL TO MINIMIZE C WORK WHEN CALLS ARE TO BE REPEATED. C W....... ARRAY CONTAINING WORKSPACE FOR COVARIANCE COMPUTATIONS. C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ C C EXTERNAL DL7VML, DV7SCP C INTEGER I, IALTM1, ICOL, ICOV, ICSET, ID, IDUM, IFAULT, IIR, 1 IOUNIT, IPCOEF, IPDUM, IPT, IPTAST, IROW, ISCALE, ISZ, 2 ITASTE, IX, J, JP, K, KP C INTEGER MAXALT, MAXAM1, LR PARAMETER (MAXALT=20, MAXAM1=MAXALT-1, LR=MAXAM1*(MAXAM1-1)/2) C DOUBLE PRECISION SCALE, SII DOUBLE PRECISION V(MAXALT), SIGMA(MAXALT,MAXALT) DOUBLE PRECISION Z(MAXAM1), SIGZ(MAXAM1,MAXAM1), R(LR) C DOUBLE PRECISION ZERO PARAMETER (ZERO=0.D0) C C SET UP V AND SIGMA MATRIX FOR MNP SPECIFICATION. C C ALTERNATIVE-SPECIFIC DUMMIES: C IALTM1 = IALT - 1 IDUM = IUSER(8) IF (IDUM.NE.0) THEN IPDUM = IUSER(12) C CASE 1: ICSET = 0. ICSET = IUSER(5) IF (ICSET.EQ.0) THEN DO 10 I = 1, IALT IX = ICDAT(I) IF (IX.NE.NALT) THEN V(I) = X(IX+IPDUM) ELSE V(I) = ZERO ENDIF 10 CONTINUE ELSE C CASE 2: ICSET.NE.0 V(IALT) = ZERO DO 20 I = 1, IALTM1 V(I) = X(I+IPDUM) 20 CONTINUE ENDIF ELSE CALL DV7SCP(IALT, V, ZERO) ENDIF C C BETA COEFFICIENTS: C NATTR = IUSER(7) IF (NATTR.NE.0) THEN IPCOEF = IUSER(11) ID = 0 DO 30 I = 1, IALT DO 30 K = 1, NATTR ID = ID + 1 V(I) = V(I) + X(IPCOEF+K)*RCDAT(ID) 30 CONTINUE ENDIF C C ALTERNATIVE-SPECIFIC ERRORS: C ICOV = IUSER(9) IF (ICOV.NE.0) THEN ICSET = IUSER(5) IF (ICSET.EQ.0) THEN DO 40 I = 1, IALT IROW = ICDAT(I) DO 40 J = 1, I ICOL = ICDAT(J) IF (ICOL.LE.IROW) THEN SIGMA(I,J) = SIGU(IROW,ICOL) ELSE SIGMA(I,J) = SIGU(ICOL,IROW) ENDIF 40 CONTINUE ELSE DO 50 I = 1, IALT DO 50 J = 1, I SIGMA(I,J) = SIGU(I,J) 50 CONTINUE ENDIF ELSE ISCALE = IUSER(18) SCALE = RUSER(ISCALE) DO 60 I = 1, IALT DO 60 J = 1, I IF (I.EQ.J) THEN SIGMA(I,J) = SCALE ELSE SIGMA(I,J) = ZERO ENDIF 60 CONTINUE ENDIF C C TASTE VARIATION: C ITASTE = IUSER(10) IF (ITASTE.EQ.1) THEN C UNCORRELATED TASTE VARIATION C SET UP W MATRIX: ID = 0 IPTAST = IUSER(14) DO 70 J = 1, IALT IPT = IPTAST DO 70 K = 1, NATTR IPT = IPT + 1 ID = ID + 1 W(K,J) = X(IPT) * RCDAT(ID) 70 CONTINUE ENDIF C IF (ITASTE.EQ.2) THEN C CORRELATED TASTE VARIATION C SET UP W MATRIX: ID = 1 IPTAST = IUSER(14) + 1 DO 80 J = 1, IALT CALL DL7VML(NATTR, W(1,J), X(IPTAST), RCDAT(ID)) ID = ID + NATTR 80 CONTINUE ENDIF IF (ITASTE.NE.0) THEN C TASTE VARIATION C ADD W(**T)W TO SIGMA: DO 100 I = 1, IALT DO 100 J = 1, I DO 90 K = 1, NATTR SIGMA(I,J) = SIGMA(I,J) + W(K,I)*W(K,J) 90 CONTINUE 100 CONTINUE ENDIF C C SYMMETRIZE SIGMA (MAY NOT BE NECESSARY???) C C IF ((ICOV.NE.0).OR.(ITASTE.NE.0)) THEN DO 110 I = 1, IALT DO 110 J = 1, I SIGMA(J,I) = SIGMA(I,J) 110 CONTINUE C ENDIF C C LOWER DIMENSION VIA STANDARD TRANSFORMATION C (REF. PAGE 43 OF DAGANZO OR BUNCH(1991)) ISZ = 0 SII = SIGMA(ICH,ICH) DO 130 JP = 1, IALT IF (JP.LT.ICH) THEN J = JP ELSE J = JP - 1 ENDIF IF (JP.NE.ICH) THEN Z(J) = V(JP)-V(ICH) DO 120 KP = 1, JP IF (KP.LT.ICH) THEN K = KP ELSE K = KP - 1 ENDIF IF(KP.NE.ICH) THEN ISZ = ISZ + 1 SIGZ(J,K)=SIGMA(JP,KP)-SIGMA(ICH,KP)-SIGMA(ICH,JP)+SII ENDIF 120 CONTINUE ENDIF 130 CONTINUE C IIR = 0 DO 150 J = 1, IALTM1 IF (SIGZ(J,J).LE.ZERO) THEN IERR = 0 RETURN ENDIF SIGZ(J,J) = SQRT(SIGZ(J,J)) Z(J) = Z(J)/SIGZ(J,J) DO 140 K = 1, J-1 IIR = IIR + 1 R(IIR) = SIGZ(J,K)/SIGZ(J,J)/SIGZ(K,K) 140 CONTINUE 150 CONTINUE C IERR = 1 CALL MNPCDF(IALTM1, Z, R, PROB, IFAULT) IF (IFAULT.NE.0) then IERR = 0 IOUNIT = IUSER(3) WRITE(IOUNIT,*) ' Problem evaluating mnpcdf' ENDIF C *** LAST LINE OF CALCP1 FOLLOWS *** END SUBROUTINE CALCDP(NPAR, X, IERR, ICH, IALT, II, ICDAT, IR, RCDAT, 1 PROB0, DP, IUSER, RUSER, MNPCDF) C C *** THIS SUBROUTINE CALCULATES FINITE-DIFFERENCE DERIVATIVES FOR *** C *** CHOICE PROBABILITIES. THIS VERSION ASSUMES THAT THE CALCPR *** C *** BEING CALLED IS THE ONE FOR MULTINOMIAL PROBIT. HOWEVER, *** C *** THE CHANGES REQUIRED FOR OTHER MODELS SHOULD BE MINOR. *** C *** NOTE: THIS SUBROUTINE REQUIRES DS7GRD, AND THE ARRAYS ALPHA *** C *** AND D SHOULD HAVE THE SAME DIMENSION AS X. *** C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ C INTEGER NPAR, IERR, ICH, IALT, II, ICDAT(*), IR, IUSER(*) DOUBLE PRECISION X(NPAR), RCDAT(*), PROB0, DP(NPAR), RUSER(*) EXTERNAL MNPCDF C EXTERNAL CALCPR, DS7GRD, DV7SCP INTEGER I, ICOV, IETA0, IPCOV, IPP, IPU, IPUP, IRC, ISCALE, ISIGP, 1 ISIGU, J, NALT, NALTM1, NFC DOUBLE PRECISION ETA, ETA0, PROB, SCALE, XTEMP C INTEGER LX DOUBLE PRECISION ONE, ZERO PARAMETER (ZERO=0.D0, ONE=1.D0, LX=60) C DOUBLE PRECISION ALPHA(LX), D(LX), WRK(6) C C *** PARAMETER USAGE *** C C SEE CALCPR AND CALCP1 C C *** BODY *** C IERR = 1 ICOV = IUSER(9) NALT = IUSER(6) NALTM1 = NALT - 1 ISCALE = IUSER(18) SCALE = RUSER(ISCALE) IETA0 = IUSER(17) ETA0 = RUSER(IETA0) C DO 10 I = 1, NPAR ALPHA(I) = ONE D(I) = ONE 10 CONTINUE C ETA = ETA0 C ETA = ETA0/PROB IRC = 0 C PROB = PROB0 20 CONTINUE CALL DS7GRD(ALPHA, D, ETA, PROB, DP, IRC, 1 NPAR, WRK, X) IF (IRC.EQ.0) GO TO 40 C IF ICOV.NE.0, SET UP AN UNPACKED SIGMA MATRIX IF (ICOV.NE.0) THEN C SQUARE THE CHOLESKY FACTOR TO GET (PACKED) SIGMA: IPCOV = IUSER(13) XTEMP = X(IPCOV) X(IPCOV) = SCALE ISIGP = IUSER(15) CALL DL7SQR(NALTM1, RUSER(ISIGP), X(IPCOV)) X(IPCOV) = XTEMP C "UNPACK" FOR EASIER ACCESS IN CALCPR: IPP = ISIGP - 1 ISIGU = IUSER(16) CALL DV7SCP(NALT*NALT, RUSER(ISIGU), ZERO) IPUP = ISIGU - 1 DO 30 I = 1, NALTM1 IPU = I + IPUP DO 30 J = 1, I IPP = IPP + 1 RUSER(IPU) = RUSER(IPP) IPU = IPU + NALT 30 CONTINUE ENDIF CALL CALCPR(NPAR, X, NFC, ICH, IALT, II, ICDAT, 1 IR, RCDAT, PROB, IUSER, RUSER, MNPCDF) IF (NFC.EQ.0) THEN IERR = 0 RETURN ENDIF GO TO 20 40 CONTINUE C C *** LAST LINE OF CALCDP FOLLOWS *** END SUBROUTINE DS7GRD (ALPHA, D, ETA0, FX, G, IRC, N, W, X) C C *** COMPUTE FINITE DIFFERENCE GRADIENT BY STWEART*S SCHEME *** C *** THIS IS SGRAD2 FROM TOMS ALGORITHM 611. C C *** PARAMETERS *** C INTEGER IRC, N DOUBLE PRECISION ALPHA(N), D(N), ETA0, FX, G(N), W(6), X(N) C C....................................................................... C C *** PURPOSE *** C C THIS SUBROUTINE USES AN EMBELLISHED FORM OF THE FINITE-DIFFER- C ENCE SCHEME PROPOSED BY STEWART (REF. 1) TO APPROXIMATE THE C GRADIENT OF THE FUNCTION F(X), WHOSE VALUES ARE SUPPLIED BY C REVERSE COMMUNICATION. C C *** PARAMETER DESCRIPTION *** C C ALPHA IN (APPROXIMATE) DIAGONAL ELEMENTS OF THE HESSIAN OF F(X). C D IN SCALE VECTOR SUCH THAT D(I)*X(I), I = 1,...,N, ARE IN C COMPARABLE UNITS. C ETA0 IN ESTIMATED BOUND ON RELATIVE ERROR IN THE FUNCTION VALUE... C (TRUE VALUE) = (COMPUTED VALUE)*(1+E), WHERE C ABS(E) .LE. ETA0. C FX I/O ON INPUT, FX MUST BE THE COMPUTED VALUE OF F(X). ON C OUTPUT WITH IRC = 0, FX HAS BEEN RESTORED TO ITS ORIGINAL C VALUE, THE ONE IT HAD WHEN DS7GRD WAS LAST CALLED WITH C IRC = 0. C G I/O ON INPUT WITH IRC = 0, G SHOULD CONTAIN AN APPROXIMATION C TO THE GRADIENT OF F NEAR X, E.G., THE GRADIENT AT THE C PREVIOUS ITERATE. WHEN DS7GRD RETURNS WITH IRC = 0, G IS C THE DESIRED FINITE-DIFFERENCE APPROXIMATION TO THE C GRADIENT AT X. C IRC I/O INPUT/RETURN CODE... BEFORE THE VERY FIRST CALL ON DS7GRD, C THE CALLER MUST SET IRC TO 0. WHENEVER DS7GRD RETURNS A C NONZERO VALUE FOR IRC, IT HAS PERTURBED SOME COMPONENT OF C X... THE CALLER SHOULD EVALUATE F(X) AND CALL DS7GRD C AGAIN WITH FX = F(X). C N IN THE NUMBER OF VARIABLES (COMPONENTS OF X) ON WHICH F C DEPENDS. C X I/O ON INPUT WITH IRC = 0, X IS THE POINT AT WHICH THE C GRADIENT OF F IS DESIRED. ON OUTPUT WITH IRC NONZERO, X C IS THE POINT AT WHICH F SHOULD BE EVALUATED. ON OUTPUT C WITH IRC = 0, X HAS BEEN RESTORED TO ITS ORIGINAL VALUE C (THE ONE IT HAD WHEN DS7GRD WAS LAST CALLED WITH IRC = 0) C AND G CONTAINS THE DESIRED GRADIENT APPROXIMATION. C W I/O WORK VECTOR OF LENGTH 6 IN WHICH DS7GRD SAVES CERTAIN C QUANTITIES WHILE THE CALLER IS EVALUATING F(X) AT A C PERTURBED X. C C *** APPLICATION AND USAGE RESTRICTIONS *** C C THIS ROUTINE IS INTENDED FOR USE WITH QUASI-NEWTON ROUTINES C FOR UNCONSTRAINED MINIMIZATION (IN WHICH CASE ALPHA COMES FROM C THE DIAGONAL OF THE QUASI-NEWTON HESSIAN APPROXIMATION). C C *** ALGORITHM NOTES *** C C THIS CODE DEPARTS FROM THE SCHEME PROPOSED BY STEWART (REF. 1) C IN ITS GUARDING AGAINST OVERLY LARGE OR SMALL STEP SIZES AND ITS C HANDLING OF SPECIAL CASES (SUCH AS ZERO COMPONENTS OF ALPHA OR G). C C *** REFERENCES *** C C 1. STEWART, G.W. (1967), A MODIFICATION OF DAVIDON*S MINIMIZATION C METHOD TO ACCEPT DIFFERENCE APPROXIMATIONS OF DERIVATIVES, C J. ASSOC. COMPUT. MACH. 14, PP. 72-83. C C *** HISTORY *** C C DESIGNED AND CODED BY DAVID M. GAY (SUMMER 1977/SUMMER 1980). C C *** GENERAL *** C C THIS ROUTINE WAS PREPARED IN CONNECTION WITH WORK SUPPORTED BY C THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS76-00324 AND C MCS-7906671. C C....................................................................... C C ***** EXTERNAL FUNCTION ***** C EXTERNAL DR7MDC DOUBLE PRECISION DR7MDC C DR7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS. C C ***** LOCAL VARIABLES ***** C INTEGER FH, FX0, HSAVE, I, XISAVE DOUBLE PRECISION AAI, AFX, AFXETA, AGI, ALPHAI, AXI, AXIBAR, 1 DISCON, ETA, GI, H, HMIN DOUBLE PRECISION C2000, FOUR, HMAX0, HMIN0, H0, MACHEP, ONE, P002, 1 THREE, TWO, ZERO C PARAMETER (C2000=2.0D+3, FOUR=4.0D+0, HMAX0=0.02D+0, HMIN0=5.0D+1, 1 ONE=1.0D+0, P002=0.002D+0, THREE=3.0D+0, 2 TWO=2.0D+0, ZERO=0.0D+0) PARAMETER (FH=3, FX0=4, HSAVE=5, XISAVE=6) C C--------------------------------- BODY ------------------------------ C IF (IRC) 50, 10, 110 C C *** FRESH START -- GET MACHINE-DEPENDENT CONSTANTS *** C C STORE MACHEP IN W(1) AND H0 IN W(2), WHERE MACHEP IS THE UNIT C ROUNDOFF (THE SMALLEST POSITIVE NUMBER SUCH THAT C 1 + MACHEP .GT. 1 AND 1 - MACHEP .LT. 1), AND H0 IS THE C SQUARE-ROOT OF MACHEP. C 10 W(1) = DR7MDC(3) W(2) = SQRT(W(1)) C W(FX0) = FX C C *** INCREMENT I AND START COMPUTING G(I) *** C 20 I = ABS(IRC) + 1 IF (I .GT. N) GO TO 120 IRC = I AFX = ABS(W(FX0)) MACHEP = W(1) H0 = W(2) HMIN = HMIN0 * MACHEP W(XISAVE) = X(I) AXI = ABS(X(I)) AXIBAR = MAX(AXI, ONE/D(I)) GI = G(I) AGI = ABS(GI) ETA = ABS(ETA0) IF (AFX .GT. ZERO) ETA = MAX(ETA, AGI*AXI*MACHEP/AFX) ALPHAI = ALPHA(I) IF (ALPHAI .EQ. ZERO) GO TO 80 IF (GI .EQ. ZERO .OR. FX .EQ. ZERO) GO TO 90 AFXETA = AFX*ETA AAI = ABS(ALPHAI) C C *** COMPUTE H = STEWART*S FORWARD-DIFFERENCE STEP SIZE. C IF (GI**2 .LE. AFXETA*AAI) GO TO 30 H = TWO* SQRT(AFXETA/AAI) H = H*(ONE - AAI*H/(THREE*AAI*H + FOUR*AGI)) GO TO 40 30 H = TWO*(AFXETA*AGI/(AAI**2))**(ONE/THREE) H = H*(ONE - TWO*AGI/(THREE*AAI*H + FOUR*AGI)) C C *** ENSURE THAT H IS NOT INSIGNIFICANTLY SMALL *** C 40 H = MAX(H, HMIN*AXIBAR) C C *** USE FORWARD DIFFERENCE IF BOUND ON TRUNCATION ERROR IS AT C *** MOST 10**-3. C IF (AAI*H .LE. P002*AGI) GO TO 70 C C *** COMPUTE H = STEWART*S STEP FOR CENTRAL DIFFERENCE. C DISCON = C2000*AFXETA H = DISCON/(AGI + SQRT(GI**2 + AAI*DISCON)) C C *** ENSURE THAT H IS NEITHER TOO SMALL NOR TOO BIG *** C H = MAX(H, HMIN*AXIBAR) IF (H .GE. HMAX0*AXIBAR) H = AXIBAR * H0**(TWO/THREE) C C *** COMPUTE CENTRAL DIFFERENCE *** C IRC = -I GO TO 100 C 50 H = -W(HSAVE) I = ABS(IRC) IF (H .GT. ZERO) GO TO 60 W(FH) = FX GO TO 100 C 60 G(I) = (W(FH) - FX) / (TWO * H) X(I) = W(XISAVE) GO TO 20 C C *** COMPUTE FORWARD DIFFERENCES IN VARIOUS CASES *** C 70 IF (H .GE. HMAX0*AXIBAR) H = H0 * AXIBAR IF (ALPHAI*GI .LT. ZERO) H = -H GO TO 100 80 H = AXIBAR GO TO 100 90 H = H0 * AXIBAR C 100 X(I) = W(XISAVE) + H W(HSAVE) = H GO TO 999 C C *** COMPUTE ACTUAL FORWARD DIFFERENCE *** C 110 G(IRC) = (FX - W(FX0)) / W(HSAVE) X(IRC) = W(XISAVE) GO TO 20 C C *** RESTORE FX AND INDICATE THAT G HAS BEEN COMPUTED *** C 120 FX = W(FX0) IRC = 0 C 999 RETURN C *** LAST CARD OF DS7GRD FOLLOWS *** END SUBROUTINE PCMRJ(NOBS, NPAR, X, NF, NEED, R, RP, UI, UR, UF) INTEGER NOBS, NPAR, NF, NEED(2), UI(*) DOUBLE PRECISION X(NPAR), R(NOBS), RP(NPAR,NOBS), UR(*) EXTERNAL UF C EXTERNAL PCMRJ1 C INTEGER ICP, IICDAT, IICH, IIIV, IIRV, IIU, INALT, IRCDAT, IRU C C *** BODY *** C IIU = UI(1) IICH = UI(2) INALT = UI(3) IIIV = UI(4) IIRV = UI(5) IICDAT = UI(6) C IRU = UI(7) ICP = UI(8) C IRW = UI(9) IRCDAT = UI(10) C CALL PCMRJ1(NOBS, NPAR, X, NF, NEED, R, RP, 1 UI(IIU), UI(IICH), UI(INALT), UI(IIIV), UI(IIRV), UI(IICDAT), 2 UR(IRU), UR(ICP), UR(IRCDAT), UF) 999 RETURN C *** LAST LINE OF PCMRJ FOLLOWS *** END SUBROUTINE PCMRJ1(NOBS, NPAR, X, NF, NEED, R, RP, 1 IUSER, ICHV, NALTV, IIV, IRV, ICDAT, 2 RUSER, CPROB, RCDAT, UF) C C *** THIS SUBROUTINE EXPANDS THE STORAGE IN UI AND UR TO MAKE THEM *** C *** COMPATIBLE WITH ESTIMATION OF CHOICE MODELS. *** C INTEGER NOBS, NPAR, NF, NEED(2), IUSER(*), ICHV(NOBS), 1 NALTV(NOBS), IIV(NOBS), IRV(NOBS), ICDAT(*) DOUBLE PRECISION X(NPAR), R(NOBS), RP(NPAR,NOBS), RUSER(*), 1 CPROB(NOBS,2), RCDAT(*) EXTERNAL UF C EXTERNAL CALCDP, CALCPR, DL7SQR, DV7SCP C INTEGER I, IALT, ICH, ICOV, IERR, II, III, IIR, IOBS, IOUNIT, 1 IPCOV, IPP, IPU, IPUP, IR, ISCALE, ISIGP, ISIGU, J, KS, 2 NALT, NALTM1, NFC DOUBLE PRECISION PROB, SCALE, XTEMP C INTEGER LX DOUBLE PRECISION ONE, ZERO PARAMETER (ZERO=0.D0, ONE=1.D0, LX=60) C C ARRAYS: C C CPROB... VECTOR FOR STORING CHOICE PROBABILITIES. CPROB(IOBS,J) C FOR J=1,2 STORES CHOICE PROBABILITIES FOR OBSERVATION C IOBS. ONE IS THE CURRENT PROBABILITY, WHILE THE OTHER C ONE IS THE PROBABILITY COMPUTED AT THE PREVIOUS TRIAL C X. THE CODE KEEPS TRACK OF WHICH IS WHICH USING THE C POINTERS STORED IN IUSER(1) AND IUSER(2). THIS IS USED C IN VARIOUS WAYS TO MAKE COMPUTATION MORE EFFICIENT. C ICHV.... VECTOR OF LENGTH NOBS. ICHV(IOBS) CONTAINS THE INDEX OF C THE CHOSEN ALTERNATIVE FOR OBSERVATION IOBS. C IIV..... VECTOR OF LENGHT NOBS. IIV(IOBS) INDICATES THE NUMBER OF C INTEGER DATA VALUES STORED IN ICDAT FOR OBSERVATION IOBS. C IRV..... VECTOR OF LENGHT NOBS. IRV(IOBS) INDICATES THE NUMBER OF C REAL DATA VALUES STORED IN RCDAT FOR OBSERVATION IOBS. C NALTV... VECTOR OF LENGHT NOBS. NALTV(IOBS) INDICATES THE NUMBER OF C CHOICES AVAILABLE FOR OBSERVATION IOBS. C C *** BODY *** C ICOV = IUSER(9) NALT = IUSER(6) NALTM1 = NALT - 1 ISCALE = IUSER(18) SCALE = RUSER(ISCALE) C IF (NEED(1).EQ.1) THEN C C *** CALCULATE RESIDUAL VECTOR *** KS = 1 IF (NEED(2).EQ.IUSER(1)) KS = 2 IUSER(KS) = NF C C IF ICOV.NE.0, SET UP AN UNPACKED SIGMA MATRIX IF (ICOV.NE.0) THEN C SQUARE THE CHOLESKY FACTOR TO GET (PACKED) SIGMA: IPCOV = IUSER(13) XTEMP = X(IPCOV) X(IPCOV) = SCALE ISIGP = IUSER(15) CALL DL7SQR(NALTM1, RUSER(ISIGP), X(IPCOV)) X(IPCOV) = XTEMP C "UNPACK" FOR EASIER ACCESS IN CALCPR: IPP = ISIGP - 1 ISIGU = IUSER(16) CALL DV7SCP(NALT*NALT, RUSER(ISIGU), ZERO) IPUP = ISIGU - 1 DO 10 I = 1, NALTM1 IPU = I + IPUP DO 10 J = 1, I IPP = IPP + 1 RUSER(IPU) = RUSER(IPP) IPU = IPU + NALT 10 CONTINUE ENDIF III = 1 IIR = 1 DO 20 IOBS = 1, NOBS ICH = ICHV(IOBS) IALT = NALTV(IOBS) II = IIV(IOBS) IR = IRV(IOBS) CALL CALCPR(NPAR, X, NFC, ICH, IALT, II, ICDAT(III), 1 IR, RCDAT(IIR), PROB, IUSER, RUSER, UF) IF ((PROB.LE.ZERO).OR.(PROB.GT.ONE).OR.(NFC.EQ.0)) THEN NF = 0 RETURN ENDIF R(IOBS) = PROB CPROB(IOBS,KS) = PROB III = III + II IIR = IIR + IR 20 CONTINUE ELSE C C *** CALCULATE JACOBIAN OF RESIDUAL VECTOR *** C KS = 1 IF (IUSER(1).NE.NF) KS = 2 IF (IUSER(KS).NE.NF) THEN IOUNIT = IUSER(3) WRITE(IOUNIT,*) ' PROBLEM WITH INITIAL ESTIMATE...' ENDIF C III = 1 IIR = 1 DO 30 IOBS = 1, NOBS ICH = ICHV(IOBS) IALT = NALTV(IOBS) II = IIV(IOBS) IR = IRV(IOBS) PROB = CPROB(IOBS,KS) CALL CALCDP(NPAR, X, IERR, ICH, IALT, II, ICDAT(III), 1 IR, RCDAT(IIR), PROB, RP(1,IOBS), IUSER, RUSER, UF) IF (IERR.EQ.0) THEN NF = 0 RETURN ENDIF III = III + II IIR = IIR + IR 30 CONTINUE ENDIF 999 RETURN C *** LAST LINE OF PCMRJ1 FOLLOWS *** END SUBROUTINE PCMRHO(NEED, F, NOBS, NF, XN, R, RD, UI, UR, W) INTEGER NEED(2), NOBS, NF, UI(*) DOUBLE PRECISION F, XN(*), R(*), RD(NOBS,*), UR(*), W(NOBS) C INTEGER ICP, IOBS, IOUNIT, IRW, WEIGHT, KS DOUBLE PRECISION OOR, VT C DOUBLE PRECISION NEGONE, ZERO PARAMETER (NEGONE=-1.D0, ZERO=0.D0) C C *** BODY *** C WEIGHT = UI(14) IF (NEED(1).EQ.1) THEN VT = ZERO IF (WEIGHT.EQ.0) THEN DO 10 IOBS = 1, NOBS VT = VT - LOG(R(IOBS)) 10 CONTINUE ELSE IRW = UI(9) DO 20 IOBS = 1, NOBS VT = VT - UR(IRW) * LOG(R(IOBS)) IRW = IRW + 1 20 CONTINUE ENDIF F = VT ELSE KS = 1 IF (UI(11).NE.NF) KS = 2 IF (UI(10+KS).NE.NF) THEN IOUNIT = UI(13) WRITE(IOUNIT,*) ' PROBLEM WITH INITIAL POINT...' NF = 0 RETURN ENDIF ICP = UI(8) IF (KS.EQ.2) ICP = ICP + NOBS IF (WEIGHT.EQ.0) THEN DO 30 IOBS = 1, NOBS OOR = NEGONE/UR(ICP) R(IOBS) = OOR W(IOBS) = R(IOBS) * OOR RD(IOBS,1) = W(IOBS) ICP = ICP + 1 30 CONTINUE ELSE IRW = UI(9) DO 40 IOBS = 1, NOBS OOR = NEGONE/UR(ICP) R(IOBS) = UR(IRW) * OOR W(IOBS) = R(IOBS) * OOR RD(IOBS,1) = W(IOBS) ICP = ICP + 1 IRW = IRW + 1 40 CONTINUE ENDIF ENDIF 999 RETURN C *** LAST LINE OF PCMRHO FOLLOWS *** END SUBROUTINE FPRINT(NOBS, NPAR, X, NF, UI, UR, UF) INTEGER NOBS, NPAR, NF, UI(*) DOUBLE PRECISION X(NPAR), UR(*) EXTERNAL UF C EXTERNAL FPRNT1 C INTEGER ICP, IICDAT, IICH, IIIV, IIRV, IIU, INALT, IRCDAT, IRU, 1 IRW C C *** BODY *** C IIU = UI(1) IICH = UI(2) INALT = UI(3) IIIV = UI(4) IIRV = UI(5) IICDAT = UI(6) C IRU = UI(7) ICP = UI(8) IRW = UI(9) IRCDAT = UI(10) C CALL FPRNT1(NOBS, NPAR, X, NF, 1 UI(IIU), UI(IICH), UI(INALT), UI(IIIV), UI(IIRV), UI(IICDAT), 2 UR(IRU), UR(IRCDAT), UR(IRW), UF) 999 RETURN C *** LAST LINE OF FPRINT FOLLOWS *** END SUBROUTINE FPRNT1(NOBS, NPAR, X, NF, 1 IUSER, ICHV, NALTV, IIV, IRV, ICDAT, 2 RUSER, RCDAT, WT, UF) C C *** THIS SUBROUTINE EXPANDS THE STORAGE IN UI AND UR TO MAKE THEM *** C *** COMPATIBLE WITH ESTIMATION OF CHOICE MODELS. *** C *** SEE PCMRJ1 DOCUMENTATION ON ARRAYS. *** C INTEGER NOBS, NPAR, NF, IUSER(*), ICHV(NOBS), NALTV(NOBS), 1 IIV(NOBS), IRV(NOBS), ICDAT(*) DOUBLE PRECISION X(NPAR), RUSER(*), RCDAT(*), WT(NOBS) EXTERNAL UF C EXTERNAL CALCPR, DL7SQR, DV7SCP C INTEGER I, IALT, ICH, ICOV, ICSET, II, III, IIR, IOBS, IOUNIT, 1 IPCOV, IPP, IPU, IPUP, IR, ISCALE, ISIGP, ISIGU, J, NALT, 2 NALTM1, NFC DOUBLE PRECISION FPROB(20), PROB, SCALE, XTEMP C INTEGER LX DOUBLE PRECISION ONE, ZERO PARAMETER (ZERO=0.D0, ONE=1.D0, LX=60) C C *** BODY *** C ICOV = IUSER(9) ICSET = IUSER(5) IOUNIT = IUSER(3) NALT = IUSER(6) NALTM1 = NALT - 1 ISCALE = IUSER(18) SCALE = RUSER(ISCALE) C WRITE(IOUNIT, 10) 10 FORMAT(//,' FINAL CHOICE SET PROBABILITIES: ',/) C C IF ICOV.NE.0, SET UP AN UNPACKED SIGMA MATRIX IF (ICOV.NE.0) THEN C SQUARE THE CHOLESKY FACTOR TO GET (PACKED) SIGMA: IPCOV = IUSER(13) XTEMP = X(IPCOV) X(IPCOV) = SCALE ISIGP = IUSER(15) CALL DL7SQR(NALTM1, RUSER(ISIGP), X(IPCOV)) X(IPCOV) = XTEMP C "UNPACK" FOR EASIER ACCESS IN CALCPR: IPP = ISIGP - 1 ISIGU = IUSER(16) CALL DV7SCP(NALT*NALT, RUSER(ISIGU), ZERO) IPUP = ISIGU - 1 DO 20 I = 1, NALTM1 IPU = I + IPUP DO 20 J = 1, I IPP = IPP + 1 RUSER(IPU) = RUSER(IPP) IPU = IPU + NALT 20 CONTINUE ENDIF III = 1 IIR = 1 DO 90 IOBS = 1, NOBS ICH = ICHV(IOBS) IALT = NALTV(IOBS) II = IIV(IOBS) IR = IRV(IOBS) DO 30 I = 1, IALT CALL CALCPR(NPAR, X, NFC, I, IALT, II, ICDAT(III), 1 IR, RCDAT(IIR), PROB, IUSER, RUSER, UF) FPROB(I) = PROB 30 CONTINUE WRITE(IOUNIT, 40) IOBS 40 FORMAT(/,' IOBS: ',I4) IF (ICSET.EQ.0) WRITE(IOUNIT,50) (ICDAT(I),I=1,IALT) 50 FORMAT(' CHOICE SET: ',20I3) WRITE(IOUNIT, 60) IALT, ICH, WT(IOBS) 60 FORMAT(' NO. OF ALTS: ',I2,' ICH: ',I2, 1 ' WT: ',F7.3) WRITE(IOUNIT, 70) (FPROB(I),I=1,IALT) 70 FORMAT(' PROBS: ',8F7.4,/,18X,8F7.4,/,18X,4F7.4) WRITE(IOUNIT, 80) FPROB(ICH) 80 FORMAT(' PROB(ICH): ',F7.4) III = III + II IIR = IIR + IR 90 CONTINUE C 999 RETURN C *** LAST LINE OF FPRNT1 FOLLOWS *** END