<div dir="ltr">FWIW, I ended up with this function:<div><br></div><div><div>(define (pick-a-maze maze-count)</div><div> (+ (if (zero? (random 2))</div><div> (/ maze-count 2)</div><div> 0)</div><div> (random-natural (/ maze-count 4))</div>
<div> (random-natural (/ maze-count 4))))</div></div><div><br></div><div style>which seems to be working out great.</div><div style><br></div><div style>Robby</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">
On Thu, Apr 11, 2013 at 12:53 PM, Joe Gilray <span dir="ltr"><<a href="mailto:jgilray@gmail.com" target="_blank">jgilray@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div dir="ltr">Really fun stuff. Thanks Neil.<div><br></div><div>-Joe</div></div><div class="gmail_extra"><br><br><div class="gmail_quote"><div><div class="h5">On Mon, Apr 8, 2013 at 2:10 PM, Robby Findler <span dir="ltr"><<a href="mailto:robby@eecs.northwestern.edu" target="_blank">robby@eecs.northwestern.edu</a>></span> wrote:<br>
</div></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div class="h5"><div dir="ltr">Thanks! What a great message!<span><font color="#888888"><div><br></div>
<div>Robby</div></font></span></div>
<div><div><div class="gmail_extra"><br><br><div class="gmail_quote">On Mon, Apr 8, 2013 at 2:47 PM, Neil Toronto <span dir="ltr"><<a href="mailto:neil.toronto@gmail.com" target="_blank">neil.toronto@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Moving this to the list because it seems general-interest-y.<div><br>
<br>
On 04/02/2013 11:08 AM, Robby Findler wrote:<br>
</div><div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi Neil; any advice on how to pick, at random, an natural number less<br>
than 604707155844263340879799987806<u></u>365649019991479003765974057968<u></u>835437533310870017962569786982<u></u>4,<br>
or other large integers?<br>
<br>
I wouldn't mind uniformly random (but if I get to make wishes, I'd like<br>
to avoid small numbers, large numbers, and numbers near n/2).<br>
<br>
I tried looking in the finite distribution and integer distribution<br>
sections of the math library but didn't see anything.<br>
</blockquote>
<br></div>
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.)<br>
<br>
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.<br>
<br>
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:<br>
<br>
#lang racket<br>
<br>
(require plot)<br>
<br>
(define is (build-list 20000 (λ (_) (random 100))))<br>
(define js (build-list 20000 (λ (_) (random 100))))<br>
(plot (density (map + is js)))<br>
<br>
<br>
Here's a more general function that chooses large natural numbers < n by adding some number of uniformly random naturals:<br>
<br>
(require math)<br>
<br>
(define random-natural-parts 3)<br>
<br>
(define (my-random-natural n)<br>
(define p (min n random-natural-parts))<br>
(define m (quotient n p))<br>
(define m-last (- n (* m (- p 1))))<br>
(+ (apply + (build-list (- p 1) (λ (_) (random-natural (+ m 1)))))<br>
(random-natural m-last)))<br>
<br>
(define xs<br>
(build-list 20000 (λ (_) (my-random-natural<br>
151223462345987123498759238861<u></u>3))))<br>
(plot (density xs))<br>
(printf "max xs = ~v~n" (apply max xs))<br>
(printf "min xs = ~v~n" (apply min xs))<br>
<br>
<br>
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.<br>
<br>
--<br>
<br>
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:<br>
<br>
(define (my-random-natural* n)<br>
(define x (my-random-natural n))<br>
(cond [(or (not (power-of-two? x)) ((random) . < . 1/4)) x]<br>
[else (my-random-natural* n)]))<br>
<br>
;; Plot a histogram<br>
(define-values (ys cs)<br>
(let*-values<br>
([(ys) (build-list 20000 (λ (_) (my-random-natural* 40)))]<br>
[(ys cs) (count-samples ys)])<br>
(sort-samples < ys cs)))<br>
<br>
(plot (discrete-histogram (map list ys cs)))<br>
<br>
<br>
Here, "rejection" means trying again, since you always want an answer.<br>
<br>
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:<br>
<br>
(define (floor-log2 x)<br>
(- (integer-length x) 1))<br>
<br>
(define (ceiling-log2 x)<br>
(integer-length (- x 1)))<br>
<br>
(define (nearest-pow2 x)<br>
(cond [(zero? x) 1]<br>
[else<br>
(define x0 (expt 2 (floor-log2 x)))<br>
(define x1 (expt 2 (ceiling-log2 x)))<br>
(if ((abs (- x x0)) . < . (abs (- x x1))) x0 x1)]))<br>
<br>
<br>
Then make and use a function that returns acceptance probabilities:<br>
<br>
(define pow2-prob 1/4)<br>
(define pow2-spread 4)<br>
<br>
(define (accept-prob x)<br>
(define pow2 (nearest-pow2 x))<br>
(define d (min 1 (/ (abs (- x pow2)) pow2-spread)))<br>
;; A convex combination of 1 and pow2-prob, with fraction d<br>
(+ d (* (- 1 d) pow2-prob)))<br>
<br>
(define (my-random-natural* n)<br>
(define x (my-random-natural n))<br>
(cond [((random) . < . (accept-prob x)) x]<br>
[else (my-random-natural* n)]))<br>
<br>
<br>
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.<br>
<br>
> (map accept-prob '(16 17 18 19 20 21))<br>
'(1/4 7/16 5/8 13/16 1 1)<br>
<br>
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.<span><font color="#888888"><br>
<br>
Neil ⊥<br>
<br>
</font></span></blockquote></div><br></div>
</div></div><br></div></div>____________________<br>
Racket Users list:<br>
<a href="http://lists.racket-lang.org/users" target="_blank">http://lists.racket-lang.org/users</a><br>
<br></blockquote></div><br></div>
</blockquote></div><br></div>