Histogram.hs

{- |
We have a set of symbols.
We draw from this set n times with replacement.
We then sort the symbols by their drawn frequency.
What is the expected distribution of the histogram?
-}
module Numeric.Probability.Example.Histogram where

import qualified Numeric.Probability.Distribution as Dist
import Control.Monad (replicateM, guard)

import qualified Data.NonEmpty as NonEmpty
import qualified Data.IntMap as IntMap
import qualified Data.Map as Map
import Data.Foldable (for_)
import Data.IntMap (IntMap)
import Data.Map (Map)

import Text.Printf (printf)

{- $setup
>>> import qualified Combinatorics as Comb
-}


example :: (Ord a) => NonEmpty.T [] a -> Int -> Dist.T Rational (Map Int Int)
example set n = Dist.norm $ do
   let x = NonEmpty.head set
   xs <- replicateM (n-1) $ Dist.uniform $ NonEmpty.flatten set
   return $ histogram $ Map.elems $ histogram (x:xs)

flattenHistogram :: (Num a) => Map Int a -> [a]
flattenHistogram histo =
   Map.elems $ Map.unionWith (+) histo $ Map.fromList $
   decorate 0 [1 .. fst $ Map.findMax histo]

{- |
>>> exampleList ('a'!:['b'..'k']) 4
fromFreqs [([4],720 % 1331),([2,1],540 % 1331),([1,0,1],40 % 1331),([0,2],30 % 1331),([0,0,0,1],1 % 1331)]
-}
exampleList :: (Ord a) => NonEmpty.T [] a -> Int -> Dist.T Rational [Int]
exampleList set n = fmap flattenHistogram $ example set n

histogram :: (Ord a) => [a] -> Map a Int
histogram = Map.fromListWith (+) . decorate 1

decorate :: b -> [a] -> [(a,b)]
decorate label = map (flip (,) label)


{-
3: https://oeis.org/A038207
4: https://oeis.org/A027465
5: https://oeis.org/A038231
6: https://oeis.org/A038243
7: https://oeis.org/A038255
8: https://oeis.org/A027466
9: https://oeis.org/A038279

expected $ example [1..k] n
j -> binomial n j * (k-1)^(n-j)

This also counts the symbols with zero occurrences.

We might prove this using a recurrence.
-}
expected :: Dist.T Rational (Map Int Int) -> Map Int Rational
expected =
   foldl (Map.unionWith (+)) Map.empty .
   map (\(x,p) -> fmap ((p*) . fromIntegral) x) .
   Dist.decons

visualize :: Rational -> Map Int Rational -> IO ()
visualize scale m =
   for_ [1 .. fst $ Map.findMax m] $ \n ->
      let freq = Map.findWithDefault 0 n m in
      printf "%6.1f %s\n" (fromRational freq :: Double) $
         replicate (round (scale * freq)) '*'



{-
https://oeis.org/A000041

>>> map (\n -> length $ partitions n n) [0..20]
-}
partitions :: Int -> Int -> [IntMap Int]
partitions maxBin =
   let go total multi =
         case compare multi 1 of
            LT -> guard (total == 0) >> [IntMap.empty]
            EQ -> guard (total <= maxBin) >> [IntMap.singleton multi total]
            GT ->
               concat $
               zipWith
                  (\j amount ->
                     map (IntMap.singleton multi j <>) $
                     go (total-amount) (multi-1))
                  [0..maxBin]
                  [0,multi..total]
   in \total -> go total total

{-
*Numeric.Probability.Visualize NE GP Comb> GP.plotList [] $ zipWith (*) (map fromInteger $ Comb.binomialSeq 110) (iterate (/27) ((27/28)**110::Double))
*Numeric.Probability.Visualize NE GP Comb> let xs = zipWith (*) (map fromInteger $ Comb.binomialSeq 110) (iterate (/27) ((27/28)**110::Double))
*Numeric.Probability.Visualize NE GP Comb> sum xs
1.0000000000000022


let histo n k = map (% k^(n-1)) $ zipWith (*) (Comb.binomialSeq n) (reverse $ take (fromInteger $ n+1) $ iterate (*(k-1)) 1)

>>> sum $ histo 110 28
28 % 1
>>> sum $ zipWith (*) [0..] $ histo 110 28
110 % 1
-}