E

August 13, 2010

We begin with a function f to count the random numbers required to reach a sum greater than one:

(define (f)
  (do ((i 0 (+ i 1))
       (s 0.0 (+ s (rand))))
      ((< 1 s) i)))

Now we run the simulation n times:

(define (e n)
  (do ((i 0 (+ i 1))
       (s 0.0 (+ s (f))))
      ((= i n) (/ s i))))

And the answer is:

> (e #e1e6)
2.718844

Run f enough times, and the surprising average is the transcendental number e. A good explanation of the math involved is given by John Walker. Mathematically, this puzzle is known as the uniform sum distribution.

We used rand from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/AWzS1d7U.

About these ads

Pages: 1 2

16 Responses to “E”

  1. [...] Praxis – E By Remco Niemeijer In today’s Programming Praxis exercise our task is to determine how many random numbers between 0 and 1 we [...]

  2. Remco Niemeijer said

    My Haskell solution (see http://bonsaicode.wordpress.com/2010/08/13/programming-praxis-e/ for a version with comments):

    import Control.Monad
    import System.Random
    
    f :: Float -> IO Int
    f s = if s > 1 then return 0 else fmap succ . f . (s +) =<< randomRIO (0, 1)
    
    test :: Int -> IO Float
    test n = fmap ((/ toEnum n) . toEnum . sum) $ replicateM n (f 0)
    
    main :: IO ()
    main = print =<< test 100000
    
  3. #include <stdio.h>
    #include <stdlib.h>
    
    main()
    {
        int trials = 1000000 ;
        int ttrials = 0;
        int i, j ;
        double sum ;
    
    
        for (i=0; i<trials; i++)
            for (sum=0., ttrials++; (sum += drand48()) < 1.; ttrials++) ;
        printf("%f trials on average.\n", (double) ttrials / trials) ;
    }
    
  4. A couple of improvements could be made. I left the j variable in accidently, it is unused. In C++, the scope of the sum value could be confined to the for loop in which it appears. A do { } while loop would mean that I didnt have two different increments of the ttrial variable.

  5. Here’s my solution in Clojure:

    (defn to-one []
      (loop [times 0 sum 0]
        (cond 
           (> sum 1) times
           :else (recur (inc times) (+ sum (rand))))))
    
    (defn to-one-avg
      [n]
      (/ (reduce + (repeatedly n to-one))
          n 1.0))
    
    user=> (to-one-avg 10000000)
    2.7184935
    
  6. Eric Hanchrow said

    PLT “racket”:

    #! /bin/sh
    #| Hey Emacs, this is -*-scheme-*- code!
    #$Id$
    exec mzscheme -l errortrace –require “$0″ –main — ${1+”$@”}
    |#

    #lang scheme
    ;; http://programmingpraxis.com/2010/08/13/e/

    (define (trial)
    (let loop ([number-of-numbers 0]
    [sum 0])
    (if (inexact
    (/ (apply + outcomes)
    (length outcomes)))))

    (provide main)

    /me tips his hat to Mark VandeWettering, Oregon ’86 or thereabouts

  7. Eric Hanchrow said

    Gaah. Stupid wordpad.

  8. Aidan said

    Here’s a Python implementation:

    import random
    import sys

    random.seed()

    numIterations = int(sys.argv[1])
    count = 0.0
    sum = 0.0

    for i in range(numIterations):
        sum += random.random()
        if sum > 1.0:
            count += 1.0
            sum = 0.0

    print numIterations/count

  9. slabounty said

    A naive ruby implementation

    count_sum = 0
    iterations = 1000
    (1..iterations).each do
        sum = 0
        count = 0
        while sum < 1.0 do
            sum += rand
            count += 1
        end
    
        count_sum += count
    end
    
    puts "#{count_sum.to_f / iterations }"
    
  10. Here’s a Java solution:

    package progpraxis;
    
    import java.util.ArrayList;
    import java.util.List;
    import java.util.Random;
    
    public class E {
    
        public static void main(String[] args) {
            Random generator = new Random();
            List<Integer> totalAttempts = new ArrayList<Integer>();
            for (int i = 1; i <= 100000; i++) {
                double total = generator.nextDouble();
                int attempts = 1;
                while (total <= 1.0) {
                    total += generator.nextDouble();
                    attempts++;
                }
                totalAttempts.add(attempts);
            }
    
            double average = 0d;
            for (Integer number : totalAttempts) {
                average += number;
            }
            average = average / totalAttempts.size();
            System.out.println("Average before surpassing 1: " + average);
        }
    }
    

    The output, as expected, is 2.72197.

  11. Axio said

    (* OCaml *)

    let test () =
      let acc = ref 0. and cpt = ref 0 in
      while !acc < 1. do
        incr cpt;
        acc := !acc +. (Random.float 1.);
      done;
      !cpt;;

    let mean n =
      let rec aux n =
        if n = 0
        then 0
        else test () + aux (n-1)
      in
      (float_of_int (aux n)) /. (float_of_int n);;

    print_string ((string_of_float (mean 10000))^”\n”);;

  12. Axio said

    (* Purely functional looks even better (well, as in “no references”…) *)

    let rec test acc cpt =
      if acc < 1.
      then test (acc +. (Random.float 1.)) (cpt + 1)
      else cpt

    let mean n =
      let rec aux n =
        if n = 0
        then 0
        else test 0. 0 + aux (n-1)
      in (float_of_int (aux n)) /. (float_of_int n);;

    print_string ((string_of_float (mean 10000))^”\n”);;

  13. [...] about computing Euler’s Number, e, using Clojure. It wasn’t until I saw the idea on Programming Praxis, though, that I decided to just do [...]

  14. Jonathan Barton said

    //quick and dirty javascript

    var getRandomCount = function(ctr,buf){
    while(buf < 1){
    var nbr = Math.random();
    ctr = ctr + 1;
    buf = buf + nbr;
    }
    return ctr;
    }

    var getAverageNumberOfTimes = function(iter){
    var total = 0;
    for(i=0;i<iter;i++){
    total = total + getRandomCount(0,0);
    }
    return total;
    }

    alert(getAverageNumberOfTimes(1000)/1000);

  15. Rodrigo Menezes said

    C#

    Random randNum = new Random();

    int trials = 1000000;
    long totalTrials = 0;
    double i = 0;

    for (int runs = 0; runs < trials; i = 0, runs++)
    for (totalTrials++ ; (i += randNum.NextDouble()) < 1.0; totalTrials++) ;

    Console.WriteLine((float)totalTrials / trials);
    Console.ReadLine();

  16. David said

    Here is a FORTH version, without floating point.

    A random 32 bit number is treated as binary random between 0-1. When overflow occurs (requires 33 bits), it is interpreted as exceeding 1.0.

    \ multiply with carry RNG
    \ quickie MWC-0  (no lag, 2^63 period is good enough)
    
    create seed  counter ,
    create carry 48313 ,
    
    : random  ( -- r )
       seed @ 4164903690 um*  carry @ 0 d+
       carry ! dup seed ! ;
    
    : trial ( -- n )
        0 0.    ( stack: count rng )
        BEGIN
            random 0 D+
            rot 1+ -rot
        dup 0> UNTIL 2drop ;
    
    : trials ( n -- )
        locals| n |
        0
        n 0 DO
            trial +
        LOOP
        10000 n */
         s>d <# # # # # [char] . hold #s #> type space ;
    

    And here is session with SwiftForth :-

    1000 trials 2.7150  ok
    10000 trials 2.7191  ok
    100000 trials 2.7163  ok
    713742 trials 2.7175  ok
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 630 other followers

%d bloggers like this: