SUBROUTINE GAUT(N,AMU,H,LDX,D,DINV,LDD,A,B,LA,LB,LE,NITER,X) C This routine draws from a truncated multivariate normal distribution C without rejection. The truncation is formed by linear restrictions C equal in number to the dimension of the distribution, +infinity and C -infinity allowed. The truncated distribution is C X ~ N(AMU, SIGMA), A < D*X < B C References: J. Geweke, 1991. "Efficient Simulation from the Multivariate C Normal and Student-t Distributions Subject to Linear C Constraints," in E.M. Keramidas (ed.), Computing Science and C Statistics: Proceedings of the 23rd Symposium on the C Interface, pp. 571-578. Fairfax, VA: Interface Foundation C of North America, Inc. C Inputs: C N Dimension of the multivariate normal distribution C AMU Mean of the distribution C H(NxN) Precision matrix of distribution C LDX Leading dimension of H, D, and DINV in calling program C D(NxN) Matrix of linear combinations of X for constraints C DINV(NxN) Inverse of D C A(N) Lower bounds for linear combinations C B(N) Upper bounds for linear combinations C LA(N) If .TRUE., no corresponding lower bound C LB(N) If .TRUE., no corresponding upper bound C LE(N) If .TRUE., equality constraint C NITER Number of Gibbs iterations to perform C Input/output: C X Random vector IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) LOGICAL LA(N),LB(N),LE(N) COMMON/SCRA/A1(LD,LD),A2(LD,LD),V1(LD),V2(LD) DIMENSION AMU(N),H(LDX,N),DINV(LDD,N),A(N),B(N),X(N) C Form mean of Z =Nu in V1 CALL DMRRRR(N,N,D,LDD,N,1,AMU,N,N,1,V1,N) C Form precision of Z in A1 CALL UXTAXF(N,N,DINV,LDD,H,LDX,A1,LD) C Retrieve current Z vector in V2 CALL DMRRRR(N,N,D,LDD,N,1,X,N,N,1,V2,N) CALL UVSUB(N,V2,V1,V2) C Make drawings for Z=D*X. DO 20 ITER=1,NITER DO 20 I=1,N IF(LE(I))GO TO 20 AA=0.0D0 DO 10 J=1,N 10 IF(I.NE.J)AA=AA-(A1(I,J)/A1(I,I))*V2(J) V2(I)=XNOT(AA,A1(I,I),A(I)-V1(I),B(I)-V1(I),LA(I),LB(I)) 20 CONTINUE C Transform from Z to X (IMSL-MATH p. 982). CALL DMRRRR(N,N,DINV,LDD,N,1,V2,N,N,1,V1,N) CALL UVADD(N,V1,AMU,X) RETURN END REAL*8 FUNCTION GAUTK(N,AMU,H,LDX,D,DINV,LDD,A,B,LA,LB,LE,NITER,X) C This function returns the log kernel of a truncated normal distribution. C If the constraints are satisfied then the kernel is the same as for the C untruncated distribution. If they are violated, a small negative number C is returned. C Inputs: C N Dimension of the multivariate normal distribution C AMU Mean of the distribution C H(NxN) Precision matrix of distribution C LDX Leading dimension of H, D, and DINV in calling program C D(NxN) Matrix of linear combinations of X for constraints C DINV(NxN) Inverse of D C A(N) Lower bounds for linear combinations C B(N) Upper bounds for linear combinations C LA(N) If .TRUE., no corresponding lower bound C LB(N) If .TRUE., no corresponding upper bound C LE(N) If .TRUE., equality constraint C NITER Number of Monte Carlo draws to use C X(N) Point of kernel evaluation C Output: C GAUTK Log kernel IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/A1(LD,LD),A2(LD,LD),V1(LD),V2(LD) DIMENSION AMU(N),H(LDX,N),DINV(LDD,N),A(N),B(N) LOGICAL LA(N),LB(N),LE(N) CALL DMURRV(N,N,D,LDD,N,X,1,N,V1) DO 10 I=1,N IF((.NOT.LA(I)).AND.(V1(I).LT.A(I)))GO TO 20 10 IF((.NOT.LB(I)).AND.(V1(I).GT.B(I)))GO TO 20 GAUTK=GAUK(N,AMU,H,LDX,X) RETURN 20 GAUTK=-700.0D0 RETURN END REAL*8 FUNCTION GAUTN(N,AMU,H,LDX,D,DINV,LDD,A,B,LA,LB,LE,NITER) C This function returns the log normalizing factor of a truncated normal C distribution. The normalilzing factor, when multiplied by the kernel C evaluated in GAUTK, yields the properly normalized probability density C of the truncated normal distribution. The normalizing factor is the C one for the untruncated normal distribution, divided by the probability C that the untruncated normal would satisfy the inequalities in the C truncated normal distribution. This probability is approximated using C the GHK simulation algorithm. C Inputs: C N Dimension of the multivariate normal distribution C AMU Mean of the distribution C H(NxN) Precision matrix of distribution C LDX Leading dimension of H, D, and DINV in calling program C D(NxN) Matrix of linear combinations of X for constraints C DINV(NxN) Inverse of D C A(N) Lower bounds for linear combinations C B(N) Upper bounds for linear combinations C LA(N) If .TRUE., no corresponding lower bound C LB(N) If .TRUE., no corresponding upper bound C LE(N) If .TRUE., equality constraint C NITER Number of Monte Carlo draws to use in the GHK algorithm C Output: C GAUTN Log normalizing factor IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(LD=100) COMMON/SCRA/A1(LD,LD),A2(LD,LD),V1(LD),V2(LD) DIMENSION AMU(N),H(LDX,N),DINV(LDD,N),A(N),B(N) GAUTN=GAUN(N,H,LDX) 1 -DLOG(GHK(N,AMU,H,LDX,D,DINV,LDD,A,B,LA,LB,LE,NITER,SE)) RETURN END REAL*8 FUNCTION GHK(N,AMU,H,LDX,D,DINV,LDD,A,B,LA,LB,LE,SD,TIME) C This routine approximates the log probability that a random variable C with the distribution X~N(AMU,H^-1) satisfies the inequality constraints C A