random walk over a unit cube

simant.fsx

// Consider an ant taking a random walk over corners of a cube.
// What are the average number of steps to arrive at the far corner?

// user interface is next line: you can change the sample size
let nSamples = 1_000_000

type Direction = X | Y | Z
/// walk in the given direction 
let walk (x,y,z) direction = 
  let step along = if along = 0 then 1 else 0
  match direction with
  | X -> step x,y,z
  | Y -> x,step y,z
  | Z -> x,y,step z

// test walk
let test1 = 
  [
    walk (0,0,0) X  =  (1,0,0)
    walk (0,0,0) Y  =  (0,1,0)
    walk (0,0,0) Z  =  (0,0,1)
    walk (0,0,1) Z  =  (0,0,0)
  ]

/// convert integer to direction
let toDirection n = if n = 0 then X else if n=1 then Y else Z

// test toDirection 
let test2 = 
  [
    toDirection 0 = X
    toDirection 1 = Y
    toDirection 2 = Z
  ]

/// random direction dispenser. not thread safe.
let ranDir seed =
    let r = System.Random(seed)
    fun() -> toDirection (r.Next(3))

/// walk from origin to opposite corner using random direction dispenser
let simAnt r = 
    let rec simAnt' here r =
        if here = (1,1,1) then 0 else 1+(simAnt' (walk here (r())) r )
    simAnt' (0,0,0) r

/// dispense fixed set of directions for testability
let fixdir items = 
    let mutable i = 0;
    fun() -> 
        let result = Array.item i items
        i <- i+1
        result

// test simulation using known directions
let test3 = [
    simAnt (fixdir [|X;Y;Z|]) = 3       // shortest journey 
    simAnt (fixdir [|X;Y;Z;Z;Y|]) = 3   // unused steps
    simAnt (fixdir [|X;Y;X;Z;X|]) = 5   // longer journey
    ]

nSamples |> printf "%A ants walking..."

// run many simulations in parallel 
// using independently seeded random number dispensers to solve
// the problem of the system random number generator not being thread safe

let randomSeeds = System.Random()
let simAntFromSeed = ranDir >> simAnt
let results = [|1..nSamples|] |> 
                Array.map (fun _ -> randomSeeds.Next()) |> 
                Array.Parallel.map simAntFromSeed

results |> Seq.averageBy float |> printfn "simple average %A" // average

// some statistics 
let median items = (Array.sort items).[items.Length/2]
let indexInProportion n p = int (float n * p)

type SomeStats = {Samples: int; Avg: float; Min:int; Max:int; Median:int; Percentile90: int}
let stats manyResults =
    let sortedResults = Array.sort manyResults;
    let n90 = indexInProportion sortedResults.Length 0.90
    {Samples=sortedResults.Length
     Avg=Seq.averageBy float manyResults
     Min=Seq.min sortedResults
     Max=Seq.max sortedResults
     Median=median sortedResults
     Percentile90=sortedResults.[n90]}

stats results |> printfn "%A"

// plot the probabilty distribution 
let histogram = results |> 
                Array.groupBy id |> 
                Array.map (fun (steps,occurances) -> steps,occurances.Length)               
let byDescendingFrequency = histogram |> Array.sortByDescending snd
let topFrequency = byDescendingFrequency.[0] |> snd
let total = histogram |> Array.map snd |> Array.reduce (+)

let probabilityDistribution = 
    byDescendingFrequency |> Array.map (fun (steps,i) -> steps, (float i)/float(total), (float i) / (float topFrequency))

let stars n = ("".PadRight n).Replace(' ','*')
let test4 = stars 4 = "****"
let starsInProportionTo fraction max = stars (int (fraction * float max))
let test5 = starsInProportionTo 0.5 10 = "*****"

printfn "time   p"
probabilityDistribution |> 
    Array.take 20 |> 
    Array.iter(fun (steps,p,topmost) -> printfn "%4d %f %s" steps p (starsInProportionTo topmost 60))