SUBROUTINE NMXD(N,NMIX,LDMIX,IDCODE,MXCODE,X,
     1                H0,HA0,S0NU0,ANU0,R,ALPHA,H,P,LT)
C  This routine draws from the posterior distribution for a discrete
C  mixture of normals density.
C  Inputs:
C  N             Sample size
C  NMIX          Number of normals in mixture
C  LDMIX         Leading dimension of H0
C  IDCODE        =1: Mean vector ALPHA in nondescending order
C                =2: Precision vector H in nondescending order
C                =3: Probability vector P in nondescending order
C                Otherwise, no identifying restrictions beyond prior
C                   distributions applied
C  MXCODE        = 1: Scale mixture of normals (means fixed at 0)
C                = 2: Mean mixture of normals (precisions all the same)
C                Otherwise, both means and precisions are free
C  X            Nx1 Data vector
C  H0           NMIX x NMIX prior precision for Alpha (ignored if MXCODE=1)
C  HA0          Product of prior precision and prior mean for Alpha
C                  (ignored if MXCODE=1)
C  S0NU0        Vector of constant parameters of chi-square priors for H
C                  (only first element if MXCODE=2)
C  ANU0         Vector of degrees-of-freedom of chi-square priors for H
C                  (only first element if MXCODE=2)
C  R            Vector of parameters of multivariate beta prior for P
C  Input/output:
C  ALPHA       NMIX x 1 vector of means.  If MXCODE not 1, NMXD
C                 exchanges values in Gibbs sampler.  If MXCODE=1,
C                 all elements of ALPHA1 must be set the same on input;
C                 values then left unchanged by NMXD
C  H           NMIX x 1 vector of precisons.  If MXCODE not 2, NMXD
C                 exchanges values in Gibbs sampler.  If MXCODE=2,
C                 H(1) is the common precision and only H(1) is exchanged
C  P           NMIX x 1 vector of probabilities.  NMXD exchanges values
C                 in Gibbs sampler
C  LT          N x 1 vector of integers between 1 and NMIX inclusive 
C                 indicating mixture status for each observation
C                 (Need not be input, since LT is drawn first)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(LD=10)
      COMMON/AMCMC/WT,PRI,PDAT,PAR(250),ITER,NPAR,LRECRD,LWRITE,LERR
      COMMON/NMXDBL/D(LD,LD),DINV(LD,LD),A(LD),B(LD),LA(LD),LB(LD),
     1              LE(LD)
      LOGICAL LA,LB,LE
      DIMENSION X(N),H0(LDMIX,NMIX),HA0(NMIX),S0NU0(NMIX),ANU0(NMIX),
     1          R(NMIX),ALPHA(NMIX),H(NMIX),P(NMIX),LT(N)
      DIMENSION A1(LD,LD),A2(LD,LD),V1(LD),V2(LD),V3(LD),NT(LD)
C  Draw latent classifications
  100 DO 110 J=1,NMIX
      V1(J)=P(J)*DSQRT(H(J))
      V2(J)=-.5D0*H(J)
  110 NT(J)=0
      DO 130 I=1,N
      DO 120 J=1,NMIX
  120 V3(J)=V1(J)*DEXP(V2(J)*(X(I)-ALPHA(J))**2)
      L=NDISC(NMIX,V3)
      LT(I)=L
  130 NT(L)=NT(L)+1
C  Draw mean vector
  200 IF(MXCODE.EQ.1)GO TO 300
      CALL UM0SET(NMIX,NMIX,A1,LD)
      CALL DSET(NMIX,0.00D0,V1,1)
      DO 210 I=1,N
      L=LT(I)
  210 V1(L)=V1(L)+H(L)*X(I)
      DO 220 J=1,NMIX
  220 A1(J,J)=H(J)*NT(J)
      NTRY=0
  230 NTRY=NTRY+1
      IF(IDCODE.EQ.1)GO TO 240
      CALL GAU2(NMIX,H0,LDMIX,A1,LD,HA0,V1,ALPHA,V2)
      GO TO 300
  240 CALL UMADD(NMIX,NMIX,H0,LDMIX,A1,LD,A2,LD)
      CALL UVADD(NMIX,HA0,V1,V2)
      CALL DLSLDS(NMIX,A2,LD,V2,V3)
      CALL GAUT(NMIX,V3,A2,LD,D,DINV,LD,A,B,LA,LB,LE,5,ALPHA)
C  Draw precision vector
  300 IF(MXCODE.EQ.2)GO TO 340
      CALL DCOPY(NMIX,S0NU0,1,V1,1)
      DO 310 I=1,N
      L=LT(I)
  310 V1(L)=V1(L)+(X(I)-ALPHA(L))**2
      DO 315 J=1,NMIX
  315 V2(J)=ANU0(J)+NT(J)
      IF(IDCODE.EQ.2)GO TO 330
      DO 320 J=1,NMIX
  320 H(J)=CHI(V1(J),V2(J))
      GO TO 400
  330 CALL CHIO(NMIX,V1,V2,1,H)
      GO TO 400
  340 AA=S0NU0(1)
      DO 350 I=1,N
  350 AA=AA+(X(I)-ALPHA(1))**2
      H(1)=CHI(AA,ANU0(1)+N)   
      CALL DSET(NMIX-1,H(1),H(2),1)
C  Draw probability vector
  400 DO 410 J=1,NMIX
  410 V1(J)=R(J)+NT(J)
  420 CALL BETAM(NMIX,V1,P)      
      IF(IDCODE.NE.3)RETURN
      DO 430 J=2,NMIX
  430 IF(P(J).LT.P(J-1))GO TO 420
      RETURN
      END