SUBROUTINE BET(N,R,P)
C  This routine generates one realization from a mutlivariate beta, or
C  Dirichlet, distribution.  The kernel density of this distribution is
C    [p1^(r1-1)]*[p2^(r2-1)]*...*[pn^(rn-1)].
C  Inputs:
C  N    Order of distributions
C  R    Entry i contains ri.
C  Outputs:
C  P    Entries (1,...,n) contain a realization of (p1,...,pn) from the
C       distribution.
C  Reference:
C  Johnson and Kotz (1972), Section 40.5
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION R(N),P(N)
      AA=0.0D0
      DO 10 I=1,N
      CALL DRNCHI(1,2.0D0*R(I),P(I))
   10 AA=AA+P(I)
      AA=1.0D0/AA
      DO 20 I=1,N
   20 P(I)=AA*P(I)
      RETURN
      END
      REAL*8 FUNCTION BETK(N,R,P)
C  This routine evaluates the log density kernel of a vector p with a 
C  multivariate Beta(n,r) distribution.  The density kernel is
C    [p1^(r1-1)]*[p2^(r2-1)]*...*[pn^(rn-1)].
C  Inputs:
C  N    Order of distributions
C  R    Entry i contains ri.
C  P    Entry i contains pi.
C  Outputs:
C  BETK Value of log density kernel
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION R(N),P(N)
      AA=0.0D0
      DO 10 I=1,N 
   10 AA=AA+(R(I)-1.0D0)*DLOG(P(I))
      BETK=AA
      RETURN
      END
      REAL*8 FUNCTION BETN(N,R)
C  This routine evaluates the log normalizing constant for the log
C  density kernel of a random variable computed in BETN.  This
C  normalizing constant, when multiplied by the kernel in BETK,
C  yields the p.d.f. of a random vector p ~ Beta(n,r).
C  Inputs:
C  N    Order of distributions
C  R    Entry i contains ri.
C  Outputs:
C  BETK Value of log density kernel
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION R(N)
      AA=0.0D0
      BB=0.0D0
      DO 10 I=1,N
      AA=AA+R(I)
   10 BB=BB+DLNGAM(R(I))
      BETN=DLNGAM(AA)-BB
      RETURN
      END