Screenshot instructions:
Windows
Mac
Red Hat Linux
Ubuntu
Click URL instructions:
Rightclick on ad, choose "Copy Link", then paste here →
(This may not be possible with some types of ads)
From: <sebastian.egner@ph...>  20030502 15:37:05
Attachments:
Message as HTML

Dear people at Schematics, Mike Sperber forwarded me a short discussion related to using pseudo random values with given probability distribution. In other words, given probabilities p[0], .., p[n1] summing to one, construct a pseudorandom process for values x in {0..n1} such that Pr{x = k} = p[k]. I would like to bring to your attention that there is a method to draw the next pseudorandom value in O(1), after having prepared a datastructure of size O(n). This method is known as the "method of aliases" by A. J. Walker (Electronic Letters 1974). It will be much faster than the code currently used in Schemathics which uses time O(n) for drawing the next random value. For convenience, I just wrote some code for this :) If there would be a good implementation for priority queues around, the preparation phase could be improved... Cheers, Sebastian. ; Random Source for Discrete Distribution ; ======================================= ; ; Sebastian.Egner@..., 2May2003 ; in R5RSScheme using ; SRFI23 (error), ; SRFI27 (random sources), ; SRFI42 (eager comprehensions) ; load this into Scheme48 0.57 with SRFI27: ; ,config ,load srfi27.scm ; ,open srfi27 srfi23 ; ,load ec.scm ; from srfi.schemers.org, SRFI42 ; (randomsourcemakediscretes s w) ; given a source s of random bits in the sense of SRFI27 ; and a vector w of n >= 1 nonnegative real numbers, ; this procedure constructs and returns a procedure rand ; such that (rand) returns the next integer x in {0..n1} ; from a pseudo random sequence with ; ; Pr{x = k} = w[k] / Sum(w[i] : i) ; ; for all k in {0..n1}. (define (randomsourcemakediscretes s w) (define (checkweights w) (if (not (and (vector? w) (>= (vectorlength w) 1) (forall?ec (:vector wi w) (and (number? wi) (not (negative? wi)))))) (error "bad vector of weights" w))) (define (makecutoffalias w) (define eps (expt 10 6)) (define n (vectorlength w)) (define cutoff (makevector n 1)) (define alias (makevector n 0)) (define avgw (/ (sumec (:vector wi w) wi) n)) (define b (vectoroflengthec n (:vector wi w) ( wi avgw))) (let loop ((sumabsb (sumec (:vector bi b) (abs bi)))) (if (< sumabsb eps) (list cutoff alias) (let ((imin 0) (bmin (vectorref b 0)) (imax 0) (bmax (vectorref b 0))) (doec (:vector bi (index i) b) (begin (if (< bi bmin) (begin (set! imin i) (set! bmin bi))) (if (> bi bmax) (begin (set! imax i) (set! bmax bi))))) (vectorset! cutoff imin (+ 1 (/ bmin avgw))) (vectorset! alias imin imax) (vectorset! b imin 0) (vectorset! b imax (+ bmin bmax)) (loop (+ sumabsb ( (abs bmin)) ( (abs bmax)) (+ (abs (vectorref b imax))))))))) (checkweights w) (let* ((ca (makecutoffalias w)) (cutoff (car ca)) (alias (cadr ca)) (randi (randomsourcemakeintegers s)) (randu (randomsourcemakereals s)) (n (vectorlength cutoff))) (lambda () (let ((i (randi n))) (if (< (randu) (vectorref cutoff i)) i (vectorref alias i)))))) ; example ;(define r1 ; (randomsourcemakediscretes ; defaultrandomsource ; (vectorec (: i 0 29) (* 0.1 (expt 0.9 i))))) ; Geo(0.1) ; ; now (r1) > x in {0..29} with Pr{x = i} = w[i]/Sum(w[j] : j) ; where w[j] = 0.1 0.9^j. ; ;(define (hist rand n iters) ; gather a histogram ; (let ((h (makevector n 0))) ; (doec (:range i iters) ; (:let x (rand)) ; (vectorset! h x (+ (vectorref h x) 1))) ; h)) ; ; ,time (hist r1 30 1000000) ; > a vector approximating the geometric distribution  Dr. Sebastian Egner Senior Scientist Philips Research Laboratories Prof. Holstlaan 4 (postbox WY63, kamer WY 655) 5656 AA Eindhoven The Netherlands tel: +31 40 2742069 *** new telephone number fax: +31 40 2742566 *** new fax number email: sebastian.egner@... 
From: <jensaxel@so...>  20030502 22:29:39

sebastian.egner@... wrote: > > Dear people at Schematics, > > Mike Sperber forwarded me a short discussion related to using > pseudo random values with given probability distribution. > In other words, given probabilities p[0], .., p[n1] summing to one, > construct a pseudorandom process for values x in {0..n1} > such that Pr{x = k} = p[k]. > > I would like to bring to your attention that there is a method to > draw the next pseudorandom value in O(1), after having > prepared a datastructure of size O(n). This method is known > as the "method of aliases" by A. J. Walker (Electronic Letters 1974). > It will be much faster than the code currently used in Schemathics > which uses time O(n) for drawing the next random value. > > For convenience, I just wrote some code for this :) That's service. I'll incorporate this as soon as possible. [When exam start for real]  Jens Axel Søgaard 
From: <jensaxel@so...>  20030504 13:04:21

sebastian.egner@... wrote: > Mike Sperber forwarded me a short discussion related to using > pseudo random values with given probability distribution. > In other words, given probabilities p[0], .., p[n1] summing to one, > construct a pseudorandom process for values x in {0..n1} > such that Pr{x = k} = p[k]. > > I would like to bring to your attention that there is a method to > draw the next pseudorandom value in O(1), after having > prepared a datastructure of size O(n). This method is known > as the "method of aliases" by A. J. Walker (Electronic Letters 1974). > It will be much faster than the code currently used in Schemathics > which uses time O(n) for drawing the next random value. For those that have Knuth, I have found the algorithm in volume 2, section 3.4.1 (page 120). He writes: However, there is actually a faster way to select x1,...,xk with arbitrarily given probabilities, based on an ingenious approach introduced by A. J. Walker [Electronics Letters 10, 8 (1974)]. Note the "ingenious". > For convenience, I just wrote some code for this :) > If there would be a good implementation for priority queues > around, the preparation phase could be improved... Is there special requirements for the priority queue? At first glance at the exercise in Knuth, it seems that one must be able to find both the minimum and the maximum element of the priority queue. Is the priority queue used after the initial preparation of the table?  Jens Axel Søgaard 
Sign up for the SourceForge newsletter:
No, thanks