/usr/share/gnudatalanguage/astrolib/randomgam.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 82 83 84 85 86 87 88 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;+
; NAME:
; RANDOMGAM
; PURPOSE:
; GENERATE GAMMA-DISTRIBUTED RANDOM VARIABLES.
;
; AUTHOR : BRANDON C. KELLY, STEWARD OBS., APRIL 2006
;
; INPUTS :
;
; SEED - THE SEED FOR THE RANDOM NUMBER GENERATOR, CAN BE UNDEFINED.
; ALPHA, BETA - THE SHAPE PARAMETERS FOR THE GAMMA DISTRIBUTION.
;
; OPTIONAL INPUTS :
;
; NRAND - THE NUMBER OF RANDOM NUMBERS TO DRAW
;-
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function randomgam, seed, alpha, beta, nrand
if n_params() lt 3 then begin
print, 'Syntax- X = randomgam( seed, alpha, beta[, nrand] )'
return, 0
endif
if alpha le 0 or beta le 0 then begin
print, 'ALPHA and BETA must both be greater than zero.'
return, 0
endif
if n_elements(nrand) eq 0 then nrand = 1
if alpha le 1 then begin
alpha = alpha + 1
alfshift = 1
endif else alfshift = 0
d = alpha - 1d / 3
c = 1 / sqrt(9 * d)
gamma = dblarr(nrand)
nempty = nrand
empty = lindgen(nrand)
repeat begin
x = randomn(seed, nempty)
v = 1 + c * x
bad = where(v le 0, nbad)
while nbad gt 0 do begin
x2 = randomn(seed, nbad)
x[bad] = x2
v[bad] = 1 + c * x2
bad2 = where(v[bad] le 0, nbad2)
if nbad2 gt 0 then bad = bad[bad2]
nbad = bad2
endwhile
v = v^3
unif = randomu(seed, nempty)
factor = 0.5 * x^2 + d - d * v + d * alog(v)
u = where( alog(unif) lt factor, nu, comp=empty1 )
if nu gt 0 then gamma[empty[u]] = d * v[u]
nempty = nempty - nu
if nempty ne 0 then empty = empty[empty1]
endrep until nempty eq 0
if alfshift then begin
alpha = alpha - 1
gamma = gamma * (randomu(seed, nrand))^(1d / alpha)
endif
gamma = gamma / beta
return, gamma
end
|