;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;+ ; 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