probability.hs

{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE FlexibleInstances #-}
import Data.List
import Control.Monad
import Data.Function
import Numeric 

formatFloatN floatNum numOfDecimals = showFFloat (Just numOfDecimals) floatNum ""

space = 10
matchResult rank = take 10 (take rank [1,1..] ++ [0,0..])

fact :: Int -> Integer
fact n = product [1..fromIntegral n]

-- | A smooth function
type F = Double -> Double

instance Num F where
  f + g = \x -> f x + g x
  f * g = \x -> f x * g x

choose :: Int -> Int -> Integer
choose n k = product (take k' [fromIntegral n, fromIntegral n-1..]) `div`fact k' where
  k' = min k (n - k)

getPosterior :: Int -> Int -> F
getPosterior wins losses = fromRational . go . toRational where
  go winProb =
    winProb ^ wins * (1 - winProb) ^ losses * fromIntegral (choose (wins + losses) wins)

integrationStep = 1e-5

-- numerically integrate from 0 to 1
integrate :: F -> Double
integrate f = sum [f x | x <- [0,integrationStep..1]]

getProbability games mustWin winProb | p <= 1 = p where
  p =
    sum [getPosterior wins (games - wins) winProb | wins <-[mustWin..games]]

expectedValue f pdf = integrate (f * pdf) / integrate pdf

uniform a b x
  | x < a = 0
  | x > b = 0
  | otherwise = 1 / (b - a)

showF f = "["++intercalate ", " [formatFloatN (f i) 2|i<-[0,0.05..1]]++"]"

result :: F -> Int -> Int -> Int -> Double
result priorWinProb wins losses toWinsWin  = integrate (probability * posterior) / integrate posterior where
   posterior = priorWinProb * getPosterior wins losses
   mustWin = toWinsWin - wins
   games = toWinsWin - losses - 1 + mustWin
   probability = getProbability games mustWin

resultNaive :: Int -> Int -> Int -> Double
resultNaive wins losses toWinsWin = probability where
   mustWin = toWinsWin - wins
   games = toWinsWin - losses - 1 + mustWin
   probability = getProbability games mustWin (fromIntegral wins / fromIntegral (wins + losses))

uniformProb = uniform 0 1

-- converts probability distribution of strength difference to
-- probability distribution of win probability
strProb :: F -> F
strProb strengthDiff = f where
  f prob
    | prob > 1 - integrationStep = 0 {- hopefully we have probability 0 of +inf strength -}
    | prob < integrationStep = 0 {- same for -inf -}
    | otherwise =
      let d1 = log (1 / (1 / (1 - (prob - integrationStep / 2)) - 1)) in
      let d2 = log (1 / (1 / (1 - (prob + integrationStep / 2)) - 1)) in
      (strengthDiff d1 + strengthDiff d2) / 2 * abs (d2 - d1) / integrationStep

gaussStr :: Double -> F
gaussStr wideness = \strength -> exp (- strength ^ 2 / wideness)

main = do
  let wins = 1
  let losses = 0
  let limit = 2
  mapM_ print $
   [ ("uniform", result uniformProb wins losses limit)
   , ("naive", resultNaive wins losses limit)
   , ("gauss str diff with wideness 0.1", result (strProb (gaussStr 0.1)) wins losses limit)
   , ("gauss str diff with wideness 1", result (strProb (gaussStr 1)) wins losses limit)
   , ("gauss str diff with wideness 2", result (strProb (gaussStr 2)) wins losses limit)
   , ("gauss str diff with wideness 3", result (strProb (gaussStr 3)) wins losses limit)
   , ("gauss str diff with wideness 10", result (strProb (gaussStr 10)) wins losses limit)
   , ("gauss str diff with wideness 100", result (strProb (gaussStr 100)) wins losses limit)
   , ("gauss str diff with wideness 10000", result (strProb (gaussStr 10000)) wins losses limit)
   ]