http://projecteuler.net/index.php?section=problems&id=84

モノポリーのどのマスに止まりやすいか、という問題。

マルコフ連鎖だとおもう。ただ、状態数が120ぐらいある。

遷移行列を求めるのが面倒だ、というか面倒そう。 

import Data.List
import Data.Array.IArray
import Control.Arrow
import Numeric.LinearAlgebra
-- dice = 6
dice = 4
len = 40
jail = 10
go = 
g2j = 30
cc = [2,17,33]
ch@[ch1,ch2,ch3]=[7,22,36]
c1 = 11
e3 = 24
h2 = 39
[r1,r2,r3] = [5,15,25]
[u1,u2] = [12,28]
xi (i,j) = i*3+j
ix =flip divMod 3
next (i,j) (d1,d2)
| j'==3   = (jail,)  -- 3 consective double
| i'==g2j = (jail,j')  -- go to jail
| i'/=g2j = (i',j')
where i' = mod (i+d1+d2) len
j' = if d1==d2 then j+1 else 
card ((i,(j,n)),p)
| elem j cc = zipWith pack [go,jail,j][p/16,p/16,7*p/8]
| elem j ch = zipWith pack
[go,jail,c1,e3,h2,r1,nextU j,mod (j-3) len,nextR j,j]
(replicate 8 (p/16)++[p/8,3*p/8])
| otherwise = [((i,(j,n)),p)]
where pack j' p' = ((i,(j',n)),p')
nextR j | j==ch1 = r2
| j==ch2 = r3
| j==ch3 = r1
nextU j | j==ch1 = u1
| j==ch2 = u2
| j==ch3 = u1
move ix = concatMap (card.f).group.sort.map (const ix&&&next ix) $range ((1,1),(dice,dice))
where f = head&&&(/(fromIntegral$dice^2)).genericLength
transM =  fromArray2D.accumArray (+)  ((,),(len*3-1,len*3-1))$moves::Matrix Double
where moves = concatMap (map toXi.move) .map ix $[..len*3-1]
toXi = first(xi***xi)
toProb ::Vector Double -> [Double]
toProb = f.toList
where f [] = []
f (a:b:c:ds) = (sum.sort$[a,b,c]):f ds
steady = toProb.fst$until enough next (fromList.(1:).replicate (len*3-1)$,constant  $len*3)
where next (sM,s) = (sM<>transM,sM)
enough (sM,s) = pnorm Infinity (sM-s) < 0.00000001
main = print.take 3.reverse.sort.zip steady $[..]

やはり予想どおり、めんどくさかった。

条件を勘違いしていたり、間違えていたりと大変だった。

行列計算にhmatrixを利用。

ちなみにdice=6で実行すると 

[(6.233011585917959e-2,10),(3.183021123222021e-2,24),(3.0887948067887767e-2,0)]

だいたいあっているでしょう。