;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;+
; NAME:
;   MULTINOM
; PURPOSE:
; SIMULATE MULTINOMIAL RANDOM VARIABLES
;
; AUTHOR : BRANDON C. KELLY, STEWARD OBS., APR 2006
;
; INPUTS :
;
;   N - THE NUMBER OF TRIALS
;   P - A K-ELEMENT VECTOR CONTAINING THE PROBABILITIES FOR EACH
;       CLASS.
;
; OPTIONAL INPUTS :
;
;   NRAND - THE NUMBER OF RANDOM VARIABLES TO DRAW
;   SEED - THE SEED FOR THE RANDOM NUMBER GENERATOR
;
; OUTPUT :
;   NRAND RANDOM DRAWS FROM A MULTINOMIAL DISTRIBUTION WITH PARAMETERS
;   N AND P.
;-
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function multinom, n, p, nrand, seed=seed

if n_params() lt 2 then begin
    print, 'Syntax- theta = multinom( n, p,[ nrand, seed=seed] )'
    return, 0
endif

k = n_elements(p)

bad = where(p lt 0 or p gt 1, nbad)
if nbad gt 0 then begin
    print, 'All element of p must be 0 <= p <= 1.'
    return, 0
endif

if n lt 1 then begin
    print, 'N must be at least 1.'
    return, 0
endif

if n_elements(nrand) eq 0 then nrand = 1

                                ;check if binomial
if k eq 2 then begin

    binom = randomu(seed, nrand, binomial=[n, p[0]], /double)
    multi = [[binom], [n - binom]]

    return, transpose(multi)

endif

multi = lonarr(k, nrand)

for i = 0L, nrand - 1 do begin
    
    multi[0,i] = randomu(seed, 1, binomial=[n, p[0]], /double)
    j = 1L
    nj = n - total(multi[0:j-1,i])

    while nj gt 0 do begin
        
        pj = p[j] / total(p[j:*])

        multi[j,i] = randomu(seed, 1, binomial=[nj,pj], /double)

        j = j + 1
        nj = n - total(multi[0:j-1,i])

    endwhile

endfor

return, multi
end