REAL*8 FUNCTION XNOT(AMU,H,A,B,LA,LB)
C  This function generates a univariate normal random variable subject
C  to the constraint that it be in an interval (a,b), where the 
C  endpoints may be finite or infinite.
C  References:  J. Geweke, 1991.  "Efficient Simulation from the 
C               Multivariate Normal and Student-t Distributions Subject 
C               to Linear Constraints," in E.M. Keramidas (ed.), 
C               Computing Science and Statistics: Proceedings of the 
C               23rd Symposium on the Interface, pp. 571-578.  Fairfax, 
C               VA: Interface Foundation of North America, Inc.
C  Inputs:
C  AMU     Mean of parent univariate normal distribution
C  H       Precision of parent univariate normal distribution
C  A, B    Endpoints of interval; A < B if LA = LB = .FALSE.
C  LA      .TRUE. if left endpoint is - infinity; in this
C          case A is ignored.
C  LB      .TRUE. if right endpoint is + infinity; in this
C          case B is ignored.
C  Output:
C  XNOT    Random variable
      IMPLICIT REAL*8 (A-H,O-Z)
      LOGICAL LA,LB,LFLIP
      DATA EPS,T1,T2,T3,T4/2.0D0,.375D0,2.18D0,.725D0,.45D0/
      F(X)=DEXP(-.5D0*X**2)
      H1=DSQRT(H)
      IF(LA.AND.LB)GO TO 160
      LFLIP=.FALSE.
      IF(LA.OR.LB)GO TO 100
      IF(B.LE.A)GO TO 170
C ******* Finite interval
      C1=H1*(A-AMU)
      C2=H1*(B-AMU)
      IF((C1*C2).GT.0.0D0)GO TO 30
C ++++ (A,B) includes 0
      IF((C1.GT.-T1).AND.(C2.LT.T1))GO TO 20
C -- F(A) or F(B) small: full normal with rejection
   10 X=DRNNOF()
      IF(X.LT.C1)GO TO 10
      IF(X.GT.C2)GO TO 10
      GO TO 150
C -- F(A) and F(B) large: uniform importance sampling
   20 CDEL=C2-C1
   25 X=C1+CDEL*DRNUNF()
      IF(DRNUNF().GT.F(X))GO TO 25
      GO TO 150
C ++++ (A,B) excludes 0
C -- Transform to both positive
   30 IF(C1.GT.0.0D0)GO TO 40
      C=C1
      C1=-C2
      C2=-C
      LFLIP=.TRUE.
   40 F1=F(C1)
      F2=F(C2)
      IF(F2.LT.EPS)GO TO 60
      IF((F1/F2).GT.T2)GO TO 60
C  -- F(A)/F(B) not large: uniform importance sampling
   50 CDEL=C2-C1
   55 X=C1+CDEL*DRNUNF()
      IF(DRNUNF().GT.(F(X)/F1))GO TO 55
      GO TO 140
   60 IF(C1.GT.T3)GO TO 80
C -- P(X>A) and F(A)/F(B) large: half-normal with rejection
   70 X=DABS(DRNNOF())
      IF(X.LT.C1)GO TO 70
      IF(X.GT.C2)GO TO 70
      GO TO 140
C -- P(X>A) small, F(A)/F(B) large: exponential importance 
C    sampling with rejection
   80 C=C2-C1
   90 Z=-DLOG(DRNUNF())/C1
      IF(Z.GT.C)GO TO 90
      IF(DRNUNF().GT.F(Z))GO TO 90
      X=C1+Z
      GO TO 140
C ****** Half-line interval
  100 C1=H1*(A-AMU)
C -- Transform to bound from below if A = -infinity
      IF(LB)GO TO 110
      C1=-H1*(B-AMU)
      LFLIP=.TRUE.
  110 IF(C1.GT.T4)GO TO 130
C -- A not large: full normal with rejection
  120 X=DRNNOF()
      IF(X.LT.C1)GO TO 120
      GO TO 140
C -- A small: exponential importance sampling
  130 Z=-DLOG(DRNUNF())/C1
      IF(DRNUNF().GT.F(Z))GO TO 130
      X=C1+Z
  140 IF(LFLIP)X=-X
  150 XNOT=AMU+X/H1
      RETURN
C ****** Whole interval
  160 XNOT=AMU+DRNNOF()/H1
      RETURN
C  ***** Singleton
  170 XNOT=A
      RETURN
      END
      REAL*8 FUNCTION XNOTK(AMU,H,A,B,LA,LB,X)
C  This function returns the logarithm of the kernel of a univariate
C  normal distribution truncated to the interval (a,b).  If x is in the
C  interval the kernel is -h*(x-mu)^2; otherwise it is zero.
C  Inputs:
C  AMU     Mean mu of parent univariate normal distribution
C  H       Precision h of parent univariate normal distribution
C  A, B    Endpoints of interval; A < B if LA = LB = .FALSE.
C  LA      .TRUE. if left endpoint is - infinity; in this
C          case A is ignored.
C  LB      .TRUE. if right endpoint is + infinity; in this
C          case B is ignored.
C  X       Point of evaluation of kernel
C  Output:
C  XNOTK    Log kernel
      IMPLICIT REAL*8 (A-H,O-Z)
      LOGICAL LA,LB
      AA=-700.0D0
      IF((.NOT.LA).AND.(X.LT.A))GO TO 10
      IF((.NOT.LB).AND.(X.GT.B))GO TO 10
      AA=-.5D0*H*(X-AMU)**2
   10 XNOTK=AA
      RETURN
      END
      REAL*8 FUNCTION XNOTN(AMU,H,A,B,LA,LB)
C***********************************************************************
C  In the calling routine a call to INIT must be made prior to calling *
C  this routine for the first time.                                    *
C***********************************************************************    
C  This function returns the logarithm of the normalizing factor of a 
C  univariate normal distribution truncated to the interval (a,b).  The
C  normalizing factor corresponds to the log-kernel evaluation in
C  XNOTK.
C  Inputs:
C  AMU     Mean mu of parent univariate normal distribution
C  H       Precision h of parent univariate normal distribution
C  A, B    Endpoints of interval; A < B if LA = LB = .FALSE.
C  LA      .TRUE. if left endpoint is - infinity; in this
C          case A is ignored.
C  LB      .TRUE. if right endpoint is + infinity; in this
C          case B is ignored.
C  Output:
C  XNOTN   Log normalizing factor
      IMPLICIT REAL*8 (A-H,O-Z)
      LOGICAL LA,LB
      COMMON/CONST/DLOG2,DLGPI,HLG2PI,PI,PII
      H1=DSQRT(H)
      AA=0.0D0
      BB=1.0D0
      FAC=-HLG2PI+.5D0*DLOG(H)
      IF(.NOT.LA)AA=DNORDF(H1*(A-AMU))
      IF(.NOT.LB)BB=DNORDF(H1*(B-AMU))
      XNOTN=FAC-DLOG(BB-AA)
      RETURN
      END