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.
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: