Galton
April 10, 2012
In the previous exercise we had one of the simulations from A. K. Dewdney’s Five Easy Pieces, and today we have another. You’ve all seen a Galton board, invented by the Victorian statistician Sir Francis Galton; a series of pegs is arranged in diagonal columns at the top, a marble is dropped into the center of the top row, and it falls through the rows, moving left or right at each peg, until it drops into a bin at the bottom of the board, one bin under each gap in the final row of pegs. Each individual marble can go left or right at each peg with equal probability, but the overall effect after several marbles are dropped is to form a bell-shaped curve in the bins at the bottom.
Your task is to write a program to simulate a Galton board, showing the bell-shaped distribution in a histogram. When you are finished you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
[…] today’s Programming Praxis exercise, our goal is to simulate a Galton board and display the frequencies of […]
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(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 "*")))))))))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.
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.
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));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: breakMy 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)) }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: