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.[0] |> 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

More information

Link:http://fssnip.net/ch
Posted:12 years ago
Author:Mathias Brandewinder
Tags: math , markov , simulation , probability