module Statistics.Distribution.Beta
( BetaDistribution(..)
, betaDistr
) where
import Numeric.SpecFunctions
import Numeric.MathFunctions.Constants (m_NaN)
import qualified Statistics.Distribution as D
import Data.Typeable
data BetaDistribution = BD
{ bdAlpha :: !Double
, bdBeta :: !Double
} deriving (Eq,Read,Show,Typeable)
betaDistr :: Double -> Double -> BetaDistribution
betaDistr a b
| a < 0 = error $ msg ++ "alpha must be positive. Got " ++ show a
| b < 0 = error $ msg ++ "beta must be positive. Got " ++ show b
| otherwise = BD a b
where msg = "Statistics.Distribution.Beta.betaDistr: "
instance D.Distribution BetaDistribution where
cumulative (BD a b) x
| x <= 0 = 0
| otherwise = incompleteBeta a b x
instance D.Mean BetaDistribution where
mean (BD a b) = a / (a + b)
instance D.MaybeMean BetaDistribution where
maybeMean = Just . D.mean
instance D.Variance BetaDistribution where
variance (BD a b) = a*b / (apb*apb*(apb+1))
where apb = a + b
invertMono :: (Double -> Double) -> Double -> Double -> Double -> Double
invertMono f l0 h0 b = go l0 h0 where
go l h
| h l < epsilon = m
| otherwise = case compare (f m) b of
LT -> go m h
EQ -> m
GT -> go l m
where m = l + (hl)/2
epsilon = 1e-12
instance D.ContDistr BetaDistribution where
density (BD a b) x
| a <= 0 || b <= 0 = m_NaN
| x <= 0 = 0
| x >= 1 = 0
| otherwise = exp $ (a1)*log x + (b1)*log (1x) logBeta a b
quantile d p
| p < 0 = error $ "probability must be positive. Got " ++ show p
| p > 1 = error $ "probability must be less than 1. Got " ++ show p
| otherwise = invertMono (D.cumulative d) 0 1 p
instance D.MaybeVariance BetaDistribution where
maybeVariance = Just . D.variance