Really fun stuff. Thanks Neil. -Joe
On Mon, Apr 8, 2013 at 2:10 PM, Robby Findler <ro...@eecs.northwestern.edu>wrote: > Thanks! What a great message! > > Robby > > > On Mon, Apr 8, 2013 at 2:47 PM, Neil Toronto <neil.toro...@gmail.com>wrote: > >> Moving this to the list because it seems general-interest-y. >> >> >> On 04/02/2013 11:08 AM, Robby Findler wrote: >> >>> Hi Neil; any advice on how to pick, at random, an natural number less >>> than 604707155844263340879799987806**365649019991479003765974057968** >>> 835437533310870017962569786982**4, >>> or other large integers? >>> >>> I wouldn't mind uniformly random (but if I get to make wishes, I'd like >>> to avoid small numbers, large numbers, and numbers near n/2). >>> >>> I tried looking in the finite distribution and integer distribution >>> sections of the math library but didn't see anything. >>> >> >> It's not in `math/distributions' because that module doesn't have integer >> distributions yet. (Well, it does, in the style of R, where the "integers" >> are actually flonums. They're fast, but probably not exact enough, or may >> not return random numbers big enough, for what you want them for.) >> >> The random big integer functions are sort of hidden in `math/base', as >> `random-natural', `random-integer' and `random-bits'. The last is the >> fastest; use it if you want a natural number less than some power of 2. >> >> The distributions for those are uniform, so to get your wishes, you'll >> have to be a bit tricky. Here's one way: add them. If you add some number p >> of uniformly distributed numbers, as p increases, the distribution of the >> sum approaches a bell shape, which sounds like what you want. For example, >> adding two results in a triangular distribution: >> >> #lang racket >> >> (require plot) >> >> (define is (build-list 20000 (λ (_) (random 100)))) >> (define js (build-list 20000 (λ (_) (random 100)))) >> (plot (density (map + is js))) >> >> >> Here's a more general function that chooses large natural numbers < n by >> adding some number of uniformly random naturals: >> >> (require math) >> >> (define random-natural-parts 3) >> >> (define (my-random-natural n) >> (define p (min n random-natural-parts)) >> (define m (quotient n p)) >> (define m-last (- n (* m (- p 1)))) >> (+ (apply + (build-list (- p 1) (λ (_) (random-natural (+ m 1))))) >> (random-natural m-last))) >> >> (define xs >> (build-list 20000 (λ (_) (my-random-natural >> 151223462345987123498759238861**3)))) >> (plot (density xs)) >> (printf "max xs = ~v~n" (apply max xs)) >> (printf "min xs = ~v~n" (apply min xs)) >> >> >> It should work for all n >= 1. Use `random-natural-parts' to control the >> shape: 1 for uniform, 2 for triangular, 3 for piecewise quadratic, etc., >> etc. The higher this number is, the less probable the numbers near 0 and n >> will be. >> >> -- >> >> To probabilistically avoid some kinds of numbers, you could use rejection >> sampling. If you wanted all powers of two to be 1/4 the probability they >> normally would be, you could accept them with probability 1/4 like this: >> >> (define (my-random-natural* n) >> (define x (my-random-natural n)) >> (cond [(or (not (power-of-two? x)) ((random) . < . 1/4)) x] >> [else (my-random-natural* n)])) >> >> ;; Plot a histogram >> (define-values (ys cs) >> (let*-values >> ([(ys) (build-list 20000 (λ (_) (my-random-natural* 40)))] >> [(ys cs) (count-samples ys)]) >> (sort-samples < ys cs))) >> >> (plot (discrete-histogram (map list ys cs))) >> >> >> Here, "rejection" means trying again, since you always want an answer. >> >> To avoid numbers *near* a power of two, make the acceptance probability a >> function of the distance to the nearest power of two. We'll need a "nearest >> power of two" function first: >> >> (define (floor-log2 x) >> (- (integer-length x) 1)) >> >> (define (ceiling-log2 x) >> (integer-length (- x 1))) >> >> (define (nearest-pow2 x) >> (cond [(zero? x) 1] >> [else >> (define x0 (expt 2 (floor-log2 x))) >> (define x1 (expt 2 (ceiling-log2 x))) >> (if ((abs (- x x0)) . < . (abs (- x x1))) x0 x1)])) >> >> >> Then make and use a function that returns acceptance probabilities: >> >> (define pow2-prob 1/4) >> (define pow2-spread 4) >> >> (define (accept-prob x) >> (define pow2 (nearest-pow2 x)) >> (define d (min 1 (/ (abs (- x pow2)) pow2-spread))) >> ;; A convex combination of 1 and pow2-prob, with fraction d >> (+ d (* (- 1 d) pow2-prob))) >> >> (define (my-random-natural* n) >> (define x (my-random-natural n)) >> (cond [((random) . < . (accept-prob x)) x] >> [else (my-random-natural* n)])) >> >> >> Here, `pow2-prob' is the probability a power of two is accepted, and >> `pow2-spread' is the minimum distance an x can be from a power of two to >> always be accepted. >> >> > (map accept-prob '(16 17 18 19 20 21)) >> '(1/4 7/16 5/8 13/16 1 1) >> >> The `accept-prob' function can be anything, as long as it returns >> probabilities, and doesn't return small probabilities so often that >> `my-random-natural*' spends a long time looping. >> >> Neil ⊥ >> >> > > ____________________ > Racket Users list: > http://lists.racket-lang.org/users > >
____________________ Racket Users list: http://lists.racket-lang.org/users
