(in-package :randist)
 
;; /* The Dirichlet probability distribution of order K-1 is 
 
;;      p(\theta_1,...,\theta_K) d\theta_1 ... d\theta_K = 
;;         (1/Z) \prod_i=1,K \theta_i^{alpha_i - 1} \delta(1 -\sum_i=1,K \theta_i)
 
;;    The normalization factor Z can be expressed in terms of gamma functions:
 
;;       Z = {\prod_i=1,K \Gamma(\alpha_i)} / {\Gamma( \sum_i=1,K \alpha_i)}  
 
;;    The K constants, \alpha_1,...,\alpha_K, must be positive. The K parameters, 
;;    \theta_1,...,\theta_K are nonnegative and sum to 1.
 
;;    The random variates are generated by sampling K values from gamma
;;    distributions with parameters a=\alpha_i, b=1, and renormalizing. 
;;    See A.M. Law, W.D. Kelton, Simulation Modeling and Analysis (1991).
 
;;    Gavin E. Crooks <gec@compbio.berkeley.edu> (2002)
;; */
 
;; void
;; gsl_ran_dirichlet (const gsl_rng * r, const size_t K,
;;                    const double alpha[], double theta[])
;; {
;;   size_t i;
;;   double norm = 0.0;
 
;;   for (i = 0; i < K; i++)
;;     {
;;       theta[i] = gsl_ran_gamma (r, alpha[i], 1.0);
;;     }
 
;;   for (i = 0; i < K; i++)
;;     {
;;       norm += theta[i];
;;     }
 
;;   for (i = 0; i < K; i++)
;;     {
;;       theta[i] /= norm;
;;     }
;; }
 
 
(defun random-dirichlet1 (alpha theta)
  (let ((norm 0.0)
	(k (array-dimension alpha 0)))
    (loop 
       for i from 0 to (1- k)
       for gamma = (random-gamma (coerce (aref alpha i)  'double-float) 1.0d0)
       sum gamma into N
       do (setf (aref theta i) gamma)
       finally (setf norm N))
    (loop 
       for i from 0 to (1- k)
       do (setf (aref theta i) (/ (aref theta i) norm)))))