/usr/share/gnudatalanguage/astrolib/multinom.pro is in gdl-astrolib 2018.02.16+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;+
; 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
|