 4 people like it.

# Simple Markov chains

Small example illustrating how to use Arrays and Sequences to simulate simple Markov chain and evaluate their stationary distribution.

 ``` 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31: 32: 33: 34: 35: 36: 37: 38: 39: 40: 41: 42: 43: 44: 45: 46: 47: 48: 49: 50: 51: 52: 53: 54: 55: 56: 57: 58: 59: 60: 61: 62: 63: 64: 65: 66: 67: 68: 69: 70: 71: 72: 73: 74: 75: 76: 77: 78: 79: 80: 81: 82: 83: 84: 85: 86: 87: 88: 89: 90: 91: 92: 93: 94: 95: 96: 97: 98: 99: 100: 101: 102: 103: 104: 105: 106: 107: 108: 109: 110: 111: 112: 113: 114: 115: ``` ``````open System // given a roll between 0 and 1 // and a distribution D of // probabilities to end up in each state // returns the index of the state let state (D: float[]) roll = let rec index cumul current = let cumul = cumul + D.[current] match (roll <= cumul) with | true -> current | false -> index cumul (current + 1) index 0.0 0 // given the transition matrix P // the index of the current state // and a random generator, // simulates what the next state is let nextState (P: float[][]) current (rng: Random) = let dist = P.[current] let roll = rng.NextDouble() state dist roll // given a transition matrix P // the index i of the initial state // and a random generator // produces a sequence of states visited let simulate (P: float[][]) i (rng: Random) = Seq.unfold (fun s -> Some(s, nextState P s rng)) i // Vector dot product let dot (V1: float[]) (V2: float[]) = Array.zip V1 V2 |> Array.map(fun (v1, v2) -> v1 * v2) |> Array.sum // Extracts the jth column vector of matrix M let column (M: float[][]) (j: int) = M |> Array.map (fun v -> v.[j]) // Given a row-vector S describing the probability // of each state and a transition matrix P, compute // the next state distribution let nextDist S P = P |> Array.mapi (fun j v -> column P j) |> Array.map(fun v -> dot v S) // Euclidean distance between 2 vectors let dist (V1: float[]) V2 = Array.zip V1 V2 |> Array.map(fun (v1, v2) -> (v1 - v2) * (v1 - v2)) |> Array.sum // Evaluate stationary distribution // by searching for a fixed point // under tolerance epsilon let stationary (P: float[][]) epsilon = let states = P. |> Array.length [| for s in 1 .. states -> 1.0 / (float)states |] // initial |> Seq.unfold (fun s -> Some((s, (nextDist s P)), (nextDist s P))) |> Seq.map (fun (s, s') -> (s', dist s s')) |> Seq.find (fun (s, d) -> d < epsilon) // Illustration // Our Markov chain models plane delays // 0 = early, 1 = on-time, 2 = delayed // more comments at // http://clear-lines.com/blog/post/Simple-Markov-chains-in-FSharp.aspx // Transition matrix: // each row corresponds to a state, // and contains an array of probabilities // to transition to each of the states. let P = [| [| 0.10; 0.85; 0.05 |]; [| 0.10; 0.75; 0.15 |]; [| 0.05; 0.60; 0.35 |] |] let rng = new Random() // how many delays in a sequence of 1000 flights? let flights = simulate P 1 rng let delays = Seq.take 1000 flights |> Seq.filter (fun i -> i = 2) |> Seq.length // stationary distribution let longterm = stationary P 0.0001 // impact of matrix modification // improve delays after delays let strat1 = [| [| 0.10; 0.85; 0.05 |]; [| 0.10; 0.75; 0.15 |]; [| 0.05; 0.61; 0.34 |] |] // improve delays after on-time let strat2 = [| [| 0.10; 0.85; 0.05 |]; [| 0.10; 0.76; 0.14 |]; [| 0.05; 0.60; 0.35 |] |] let impact1 = stationary strat1 0.0001 let impact2 = stationary strat2 0.0001 ``````
namespace System
val state : D:float [] -> roll:float -> int

Full name: Script.state
val D : float []
Multiple items
val float : value:'T -> float (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.float

--------------------
type float = Double

Full name: Microsoft.FSharp.Core.float

--------------------
type float<'Measure> = float

Full name: Microsoft.FSharp.Core.float<_>
val roll : float
val index : (float -> int -> int)
val cumul : float
val current : int
val nextState : P:float [] [] -> current:int -> rng:Random -> int

Full name: Script.nextState
val P : float [] []
val rng : Random
Multiple items
type Random =
new : unit -> Random + 1 overload
member Next : unit -> int + 2 overloads
member NextBytes : buffer:byte[] -> unit
member NextDouble : unit -> float

Full name: System.Random

--------------------
Random() : unit
Random(Seed: int) : unit
val dist : float []
Random.NextDouble() : float
val simulate : P:float [] [] -> i:int -> rng:Random -> seq<int>

Full name: Script.simulate
val i : int
module Seq

from Microsoft.FSharp.Collections
val unfold : generator:('State -> ('T * 'State) option) -> state:'State -> seq<'T>

Full name: Microsoft.FSharp.Collections.Seq.unfold
val s : int
union case Option.Some: Value: 'T -> Option<'T>
val dot : V1:float [] -> V2:float [] -> float

Full name: Script.dot
val V1 : float []
val V2 : float []
type Array =
member Clone : unit -> obj
member CopyTo : array:Array * index:int -> unit + 1 overload
member GetEnumerator : unit -> IEnumerator
member GetLength : dimension:int -> int
member GetLongLength : dimension:int -> int64
member GetLowerBound : dimension:int -> int
member GetUpperBound : dimension:int -> int
member GetValue : [<ParamArray>] indices:int[] -> obj + 7 overloads
member Initialize : unit -> unit
member IsFixedSize : bool
...

Full name: System.Array
val zip : array1:'T1 [] -> array2:'T2 [] -> ('T1 * 'T2) []

Full name: Microsoft.FSharp.Collections.Array.zip
val map : mapping:('T -> 'U) -> array:'T [] -> 'U []

Full name: Microsoft.FSharp.Collections.Array.map
val v1 : float
val v2 : float
val sum : array:'T [] -> 'T (requires member ( + ) and member get_Zero)

Full name: Microsoft.FSharp.Collections.Array.sum
val column : M:float [] [] -> j:int -> float []

Full name: Script.column
val M : float [] []
val j : int
Multiple items
val int : value:'T -> int (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.int

--------------------
type int = int32

Full name: Microsoft.FSharp.Core.int

--------------------
type int<'Measure> = int

Full name: Microsoft.FSharp.Core.int<_>
val v : float []
val nextDist : S:float [] -> P:float [] [] -> float []

Full name: Script.nextDist
val S : float []
val mapi : mapping:(int -> 'T -> 'U) -> array:'T [] -> 'U []

Full name: Microsoft.FSharp.Collections.Array.mapi
val dist : V1:float [] -> V2:float [] -> float

Full name: Script.dist
val stationary : P:float [] [] -> epsilon:float -> float [] * float

Full name: Script.stationary
val epsilon : float
val states : int
val length : array:'T [] -> int

Full name: Microsoft.FSharp.Collections.Array.length
val s : float []
val map : mapping:('T -> 'U) -> source:seq<'T> -> seq<'U>

Full name: Microsoft.FSharp.Collections.Seq.map
val s' : float []
val find : predicate:('T -> bool) -> source:seq<'T> -> 'T

Full name: Microsoft.FSharp.Collections.Seq.find
val d : float
val P : float [] []

Full name: Script.P
val rng : Random

Full name: Script.rng
val flights : seq<int>

Full name: Script.flights
val delays : int

Full name: Script.delays
val take : count:int -> source:seq<'T> -> seq<'T>

Full name: Microsoft.FSharp.Collections.Seq.take
val filter : predicate:('T -> bool) -> source:seq<'T> -> seq<'T>

Full name: Microsoft.FSharp.Collections.Seq.filter
val length : source:seq<'T> -> int

Full name: Microsoft.FSharp.Collections.Seq.length
val longterm : float [] * float

Full name: Script.longterm
val strat1 : float [] []

Full name: Script.strat1
val strat2 : float [] []

Full name: Script.strat2
val impact1 : float [] * float

Full name: Script.impact1
val impact2 : float [] * float

Full name: Script.impact2