Numeric.Probability.Example.Kruskal

Given a row of n (~50) dice, two players start with a random dice within the first m (~5) dice. Every player moves along the row, according to the pips on the dice. They stop if a move would exceed the row. What is the probability that they stop at the same die? (It is close to one.)

Kruskal's trick: http://faculty.uml.edu/rmontenegro/research/kruskal_count/kruskal.html

Wuerfelschlange (german): http://www.math.de/exponate/wuerfelschlange.html/

Documentation

type Die = Int

type Probability = Rational

type Dist = T Probability

die :: (C prob experiment, Fractional prob) => Score -> experiment Die

type Score = Int

game :: (C prob experiment, Fractional prob) => Score -> Score -> (Score, Score) -> experiment (Maybe Score)

We reformulate the problem to the following game: There are two players, both of them collect a number of points. In every round the player with the smaller score throws a die and adds the pips to his score. If the two players somewhen get the same score, then the game ends and the score is the result of the game (Just score). If one of the players exceeds the maximum score n, then the game stops and players lose (Nothing).

gameRound :: Score -> Score -> Dist (Either (Maybe Score) (Score, Score)) -> Dist (Either (Maybe Score) (Score, Score))

gameFast :: Score -> Score -> Dist (Score, Score) -> Dist (Maybe Score)

gameFastEither :: Score -> Score -> Dist (Score, Score) -> Dist (Maybe Score)

gameFastFix :: Score -> Score -> Dist (Score, Score) -> Dist (Maybe Score)

This version could be generalized to both Random and Distribution monad while remaining efficient.

gameLeastScore :: Score -> Score -> Dist (Score, Score) -> Dist (Maybe Score)

In gameFastFix we group the scores by rounds. This leads to a growing probability distribution, but we do not need the round number. We could process the game in a different way: We only consider the game states where the lower score matches the round number.

flattenedMatrix :: Score -> [Int]

gameLeastScore can be written in terms of a matrix power. For n pips we need a n² × n² matrix. Using symmetries, we reduce it to a square matrix with size n·(n+1)/2.

p[n+1,(n+1,n+1)] p[n,(n+0,n+0)] | p[n+1,(n+1,n+2)] | | p[n,(n+0,n+1)] | | p[n+1,(n+1,n+3)] | | p[n,(n+0,n+2)] | | ... | | ... | | p[n+1,(n+1,n+6)] | = M/6 · | p[n,(n+0,n+5)] | | p[n+1,(n+2,n+2)] | | p[n,(n+1,n+1)] | | ... | | ... | | p[n+1,(n+2,n+6)] | | p[n,(n+1,n+5)] | | ... | | ... | p[n+1,(n+6,n+6)] p[n,(n+5,n+5)]

M[(n+1,(x,y)),(n,(x,y))] = 6

M[(n+1,(min y (n+d), max y (n+d))), (n,(n,y))] = 1

M[(n+1,(x1,y1)),(n,(x0,y0))] = 0

startVector :: Score -> [Int]

compareMaybe :: Ord a => Maybe a -> Maybe a -> Ordering

cumulate :: Ord a => Dist (Maybe a) -> [(Maybe a, Probability)]

runExact :: Score -> IO ()

trace :: Score -> [Score] -> [Score]

chop :: [Score] -> [[Score]]

meeting :: [Score] -> [Score] -> Maybe Score

bruteforce :: Score -> Score -> (Score, Score) -> T (Maybe Score)

This is a bruteforce implementation of the original game: We just roll the die maxScore times and then jump from die to die according to the number of pips.

runSimulation :: Score -> IO ()

latexDie :: Score -> String

latexMarkedDie :: Score -> String

latexFromChain :: [Score] -> String

latexChoppedFromChain :: [Score] -> String

makeChains :: IO ()