;+ ; NAME: ; random_non_uni ; ; PURPOSE: ; Returns a random number with a given probability distribution ; ; CALLING SEQUENCE: ; result = random_non_uni(pop, num_wanted) ; ; INPUTS: ; pop: the population of values that gives the probability distribution ; num_wanted: the number of random values to return ; ; RESTRICTIONS: ; Does NOT return multi-dimensional arrays, like randomu and randomn. ; Also, this routine doesn't do any checking of the input values, so make ; sure you know what you're doing. ; ; EXAMPLE: ; ;Some arbitrary, random population of numbers ; pop = exp(-randomu(seed, 101)) ; res = random_non_uni(pop, 10000) ; ; MODIFICATION HISTORY: ; Written by Brian Jackson (decaelus@gmail.com), 2012 June 19 ; ;- function random_non_uni, pop, num_wanted srt_pop = pop[sort(pop)] uniq_pop = srt_pop[uniq(srt_pop)] ;Calculate cumulative histogram num = double(n_elements(uniq_pop)) cum_hist = dblarr(num) cum_hist[0] = 1./num for i = 1, n_elements(cum_hist)-1 do \$ cum_hist[i] = (n_elements(where(srt_pop le uniq_pop[i])))/double(n_elements(pop)) ;To understand how this calculation works, check out: ; _Numerical Recipes in C_ (1992), ch. 7.2. rand = randomu(seed, num_wanted) ret = dblarr(num_wanted) ind = value_locate(cum_hist, rand) for i = 0, num_wanted-1 do begin aind = double(cum_hist[ind[i]+1]-rand[i])/double(cum_hist[ind[i]+1]-cum_hist[ind[i]]) ret[i] = interpol([uniq_pop[ind[i]], uniq_pop[ind[i]+1]], [0, 1], aind) end;for i return, ret end;function