probability 0.2.8.1 → 0.2.9

raw patch · 4 files changed

+122/−7 lines, 4 filesdep +non-empty

Dependencies added: non-empty

Files

probability.cabal +4/−2
src/Numeric/Probability/Example/Dice.hs +1/−1
src/Numeric/Probability/Example/DiceAccum.hs +1/−4
src/Numeric/Probability/Example/Histogram.hs +116/−0

probability.cabal view

@@ -1,6 +1,6 @@ Cabal-Version:      1.24 Name:               probability-Version:            0.2.8.1+Version:            0.2.9 License:            BSD3 License-File:       COPYRIGHT Author:             Martin Erwig <erwig@eecs.oregonstate.edu>, Steve Kollmansberger@@ -30,7 +30,7 @@ Source-Repository this   type:     darcs   location: http://code.haskell.org/~thielema/probability/-  tag:      0.2.8.1+  tag:      0.2.9   Flag splitBase@@ -38,6 +38,7 @@  Library   Build-Depends:+    non-empty >=0.3 && <0.4,     utility-ht >=0.0.12 && <0.1,     transformers >=0.4 && <0.7   If flag(splitBase)@@ -75,6 +76,7 @@     Numeric.Probability.Example.Kruskal     Numeric.Probability.Example.MontyHall     Numeric.Probability.Example.NBoys+    Numeric.Probability.Example.Histogram     Numeric.Probability.Example.Predator     Numeric.Probability.Example.Profession     Numeric.Probability.Example.Queuing

src/Numeric/Probability/Example/Dice.hs view

@@ -10,7 +10,7 @@ type Probability = Rational type Dist = Dist.T Probability -die :: Dist Die+die :: Fractional prob => Dist.T prob Die die = Dist.uniform [1..6]  -- | product of independent distributions

src/Numeric/Probability/Example/DiceAccum.hs view

@@ -14,16 +14,13 @@ import qualified Numeric.Probability.Monad as MonadExt import Numeric.Probability.Trace (Trace) -import Numeric.Probability.Example.Dice (Die, )+import Numeric.Probability.Example.Dice (die)  import Data.Maybe.HT (toMaybe)   type Score = Int --die :: Fractional prob => Dist.T prob Die-die = Dist.uniform [1..6]  roll :: Fractional prob => Trans.T prob (Maybe Score) roll =

+ src/Numeric/Probability/Example/Histogram.hs view

@@ -0,0 +1,116 @@+{- |+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+-}