src/Numeric/Probability.hs

   1 -- For `ProbabilityScale`
   2 {-# LANGUAGE DataKinds #-}
   3
   4 module Numeric.Probability (
   5   Probability,
   6   ProbabilityScale,
   7   ProbabilityBounded (..),
   8   probability,
   9   safeProbability,
  10   runProbability,
  11   assertProbability,
  12   proba0,
  13   proba1,
  14 )
  15 where
  16
  17 import Control.Monad (Monad (..))
  18 import Data.Eq (Eq)
  19 import Data.Function (id)
  20 import Data.Maybe (Maybe (..), fromJust)
  21 import Data.Monoid (Monoid (..))
  22 import Data.Ord (Ord (..))
  23 import Data.Proxy (Proxy (..))
  24 import Data.Semigroup (Semigroup (..))
  25 import Data.Validity (Validity (..), declare)
  26 import Data.Word (Word64)
  27 import GHC.Generics (Generic)
  28 import GHC.TypeNats (Natural, natVal)
  29 import Logic
  30 import Logic.Theory.Bool (type (&&))
  31 import Logic.Theory.Ord (type (<=))
  32 import Numeric.Decimal (Decimal (..), MonadThrow (..))
  33 import Numeric.Decimal qualified as Decimal
  34 import System.Random (Random)
  35 import Text.Show (Show (show))
  36 import Prelude (Bounded (..), Enum, Fractional (..), Integral, Num (..), Rational, Real, error, (^))
  37
  38 type Probability = Decimal Decimal.RoundHalfEven ProbabilityScale ProbabilityBounded
  39 instance Validity Probability where
  40   validate (Decimal wb) = validate wb
  41
  42 probability :: MonadThrow m => Rational -> m Probability
  43 probability = Decimal.fromRationalDecimalBoundedWithRounding
  44 {-# INLINE probability #-}
  45
  46 proba0 :: Probability
  47 proba1 :: Probability
  48 proba0 = fromJust (probability 0)
  49 proba1 = fromJust (probability 1)
  50
  51 safeProbability :: r ::: Rational / r <= 0 && r <= 1 -> Probability
  52 safeProbability (Named r) = fromJust (probability r)
  53
  54 runProbability :: Probability -> Rational
  55 runProbability = Decimal.toRationalDecimal
  56 {-# INLINE runProbability #-}
  57
  58 assertProbability :: Rational -> Probability
  59 assertProbability r = case probability r of
  60   Just p -> p
  61   Nothing -> error ("assertProbability: " <> show r)
  62
  63 instance Num (Decimal.Arith Probability) where
  64   (+) = Decimal.bindM2 Decimal.plusDecimalBounded
  65   (-) = Decimal.bindM2 Decimal.minusDecimalBounded
  66   (*) = Decimal.bindM2 Decimal.timesDecimalBoundedWithRounding
  67   abs = id
  68   signum m = m >>= Decimal.signumDecimalBounded
  69   fromInteger = Decimal.fromIntegerDecimalBoundedIntegral
  70
  71 instance Fractional (Decimal.Arith Probability) where
  72   (/) = Decimal.bindM2 Decimal.divideDecimalBoundedWithRounding
  73   fromRational = probability
  74
  75 -- >>> 10^19 <= (fromIntegral (maxBound :: Word64) :: Integer
  76 -- True
  77 -- >>> 10^20 <= (fromIntegral (maxBound :: Word64) :: Integer
  78 -- False
  79 type ProbabilityScale = 19
  80
  81 newtype ProbabilityBounded = ProbabilityBounded {unProbabilityBounded :: Word64}
  82   deriving (Show, Eq, Ord, Num, Real, Integral, Enum, Random, Generic)
  83 instance Bounded ProbabilityBounded where
  84   minBound = ProbabilityBounded 0
  85   maxBound = ProbabilityBounded (10 ^ (natVal (Proxy @ProbabilityScale)))
  86 instance Validity ProbabilityBounded where
  87   validate (ProbabilityBounded w) =
  88     mconcat
  89       [ declare ("The contained word is smaller or equal to 10 ^ ProbabilityScale = " <> show (10 ^ n :: Natural)) (w <= 10 ^ n)
  90       ]
  91     where
  92       n :: Natural
  93       n = natVal (Proxy @ProbabilityScale)