Galton

April 10, 2012

The program decomposes into three pieces: drop one marble, drop a series of marbles, and draw a histogram. Here’s the function to drop one marble, which merely needs to count the number of times the marble moves right, signalled when the random integer is 1:

(define (drop bins)
  (do ((b bins (- b 1))
       (z 0 (+ z (randint 2))))
      ((= b 1) z)))

Function galton drops a series of marbles:

(define (galton marbles bins)
  (let ((bs (make-vector bins 0)))
    (do ((m 0 (+ m 1))) ((= m marbles) bs)
      (let ((b (drop bins)))
        (vector-set! bs b
          (+ (vector-ref bs b) 1))))))

Then the histogram, which scales itself and is displayed horizontally, using an auxiliary function to draw the row of asterisks:

(define (stars n)
  (do ((n n (- n 1))) ((zero? n) (newline))
    (display #\*)))

(define (histogram xs)
  (let* ((max-x (apply max (vector->list xs)))
         (len (vector-length xs))
         (n (ceiling (/ max-x len))))
    (do ((x 0 (+ x 1))) ((= x len))
      (stars (ceiling (/ (vector-ref xs x) n))))))

Here’s a sample run:

> (histogram (galton 1000 8))
*
**
******
********
********
*****
**
*

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

About these ads

Pages: 1 2

8 Responses to “Galton”

  1. […] today’s Programming Praxis exercise, our goal is to simulate a Galton board and display the frequencies of […]

  2. My Haskell solution (see http://bonsaicode.wordpress.com/2012/04/10/programming-praxis-galton/ for a version with comments):

    import Control.Applicative
    import Control.Monad
    import System.Random
    
    marble :: Int -> IO Int
    marble bins = sum . take (bins - 1) . randomRs (0, 1) <$> newStdGen
    
    marbles :: Num a => Int -> Int -> IO [a]
    marbles n bins = flip fmap (replicateM n $ marble bins) (\results ->
        map (fromIntegral . length . flip filter results . (==)) [0..bins - 1])
    
    histogram :: RealFrac a => a -> [a] -> IO ()
    histogram w cols = mapM_ (\n -> 
        putStrLn $ replicate (ceiling $ n * w / maximum cols) '*') cols
    
    main :: IO ()
    main = histogram 20 =<< marbles 1000 8
    
  3. uberlazy said
    (defn 
      directions
      []
      (for [_ (range)]
        (rand-int 2)))
    
    (defn single-pass [size]
      (reduce + 0 (take size (directions))))
    
    (defn galton [size passes]
      (apply merge-with +
             (for [_ (range passes)]
               {(single-pass size) 1})))
    
    (defn histogram
      ([size passes] (histogram size passes 20))
      ([size passes scale-factor]
         (let [g (galton size passes)]
           (doseq [i (range (inc size))
                   :let [cnt (g i 0)]]
             (println (format "%5d %5d %s"
                              i
                              cnt
                              (apply str
                                     (take  (* scale-factor (/ cnt passes))
                                            (repeat "*")))))))))
    
  4. Jussi Piitulainen said

    In R, after realizing that the sum of 200 Bernoulli trials with p=0.5 (each success means falling to the right at a pin, incrementing the bin number by 1) is a binomial trial with n=200 and p=0.5, one obtains the final result thus.

    plot(table(rbinom(1000, 200, 0.5)))
    

    That is a bar plot for 1000 balls through 200 levels of pins on the board, hence 201 bins. The real histogram command of R has also the tabulation of values built-in, but I think the bar plot obtained above looks nice for this.

    hist(rbinom(1000, 200, 0.5))
    [/source]
    It's equally easy to obtain 200 Bernoulli values and sum them to get the bin number of one ball, but to repeat that for 1000 balls would require some sort of loop, I think.
    [sourcecode lang="css"]
    sum(rbinom(200, 1, 0.5))
    
  5. ardnew said

    Perl solution, output scaled to $MAX_WIDTH characters

    use strict;
    use warnings;
    
    my $MAX_WIDTH = 320;
    my $SYMBOL    = '*';
    
    sub drop
    {
      my ($nbins, $nballs, @bins) = @_;
      map { ++$bins[grep { int(rand(2)) } 1 .. $nbins] } 1 .. $nballs;
      return @bins;
    }
    
    die "\nusage:\n\t$0 <bins> <balls>\n" unless 
      scalar @ARGV > 1;
      
    printf("%5d: %s\n", $_, $SYMBOL x ($_ * $MAX_WIDTH / $ARGV[1])) 
      foreach (drop(@ARGV));
    
  6. Mike said

    Here’s a quick and dirty Galton board animation in Python 2.7.

    from collections import defaultdict
    from math import sin, cos, radians as rads
    from random import choice
    from turtle import *
    
    R = 5
    
    def galton(nlevels=11, trace=False, speed=6):
        setup(width=500, height=800, startx=0, starty=0)
        delay(0)
        pen(speed=10, shown=False)
        goto(-200, -350)
        clear()
        goto( 200, -350)
        pen(pendown=False)
        for n,y in enumerate(range(nlevels*14, -1, -14), 1):
            for x in range(-7*n, 7*n+1, 14):
                goto(x, y)
                dot(3)
        pen(speed=speed)
        
        angles = (rads(d) for d in range(0,361,36))
        addshape('marble', tuple((R*sin(r),R*cos(r)) for r in angles))
        shape('marble')
        
        hist = defaultdict(lambda:-343)
        while True:
            pen(shown=False, pendown=False)
            x,y = 0, nlevels*14
            goto(x, y)
            pen(shown=True, pendown=trace)
            delay(5)
            
            while y > 0:
                if y%14:
                    x += choice((-7,7))
                    goto(x, y)
                y -= 7
                goto(x, y)
    
            delay(1)
            while y > hist[x]:
                y -= 7
                goto(x, y)
                
            stamp()
            hist[x] += 7
            
            if hist[x] > 0:
                break
    
  7. My answer in Go:

    package main
    
    import (
    	"fmt"
    	"time"
    	"math/rand"
    )
    
    func galtonStep(ball, peg chan int) {
    	rng := rand.New(rand.NewSource(time.Now().UnixNano()))
    	for {
    		pos := <- ball
    		switch coin := rng.Float32(); true {
    			case coin < 0.5:
    				peg <- pos
    			case coin >= 0.5:
    				peg <- pos + 1
    		}
    	}
    }
    
    func galtonBoard(depth int) (ball, bin chan int) {
    	c := make(chan int)
    	ball, bin = c, c
    
    	for i := 0; i < depth; i++ {
    		prev := bin
    		bin = make(chan int)
    		go galtonStep(prev, bin)
    	}
    
    	return
    }
    
    func galtonSimulate(depth, count int) (bins []int) {
    	bins = make([]int, depth + 1)
    
    	ball, bin := galtonBoard(depth)
    	for i := 0; i < count; i++ {
    		ball <- 0
    		bins[<- bin]++
    	}
    
    	return
    }
    
    func max(values []int) int {
    	r := values[0]
    	for _, v := range(values) {
    		if v > r {
    			r = v
    		}
    	}
    	return r
    }
    
    func printHist(cols int, hist []int) {
    	maxv := max(hist)
    
    	for bin, v := range(hist) {
    		count := cols * v / maxv
    		fmt.Printf("%02d: ", bin)
    		for i := 0; i < count; i++ {
    			fmt.Print("*")
    		}
    		fmt.Println()
    	}
    }
    
    func main() {
    	printHist(65, galtonSimulate(20, 1000))
    }
    
  8. David said

    Forth version, with a RNG developed in an earlier exercise (not included here.)

    include random
    
    9 constant #Bins
    25 constant HGRAM_WIDTH
    
    : int-array ( n -- )
        create 
            cells allot
        does>
            swap cells + ;
    
    #Bins int-array galton-board
    
    : init-galton  ( -- )  \ clear the marbles from the machine
        ['] galton-board >body  #bins cells  0 fill ;
    
    : .galton
        #bins 0 DO  i galton-board ?  LOOP ;
    
    : drop-marble  ( -- )   \ drop a marble into the machine
        0
        [ #bins 1- ] literal 0 DO  2 random +  LOOP
        galton-board ++ ;
    
    : max-galton  ( -- n )
        0
        #bins 0 DO  i galton-board @ max  LOOP ;
    
    : stars  ( n -- )
        0 ?DO  [char] * emit  LOOP ;
    
    : histogram  ( width -- )
        max-galton  locals| board-max column-max |
        cr ." BIN"
    
        #bins 0 DO
            cr
            i 3 .r space
            i galton-board @  column-max board-max */  stars
         LOOP ;
    
    : galton ( marbles )
        init-galton
        0 DO  drop-marble  LOOP
        HGRAM_WIDTH histogram ;
    

    Execution:

    100000 galton
    BIN
      0
      1 **
      2 *********
      3 ********************
      4 *************************
      5 ********************
      6 **********
      7 **
      8  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 645 other followers

%d bloggers like this: