1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
|
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module : Statistics.Distribution.ChiSquared
-- Copyright : (c) 2010 Alexey Khudyakov
-- License : BSD3
--
-- Maintainer : bos@serpentine.com
-- Stability : experimental
-- Portability : portable
--
-- The chi-squared distribution. This is a continuous probability
-- distribution of sum of squares of k independent standard normal
-- distributions. It's commonly used in statistical tests
module Statistics.Distribution.ChiSquared (
ChiSquared
, chiSquaredNDF
-- * Constructors
, chiSquared
, chiSquaredE
) where
import Control.Applicative
import Data.Aeson (FromJSON(..), ToJSON, Value(..), (.:))
import Data.Binary (Binary(..))
import Data.Data (Data, Typeable)
import GHC.Generics (Generic)
import Numeric.SpecFunctions ( incompleteGamma,invIncompleteGamma,logGamma,digamma)
import Numeric.MathFunctions.Constants (m_neg_inf)
import qualified System.Random.MWC.Distributions as MWC
import qualified Statistics.Distribution as D
import Statistics.Internal
-- | Chi-squared distribution
newtype ChiSquared = ChiSquared
{ chiSquaredNDF :: Int
-- ^ Get number of degrees of freedom
}
deriving (Eq, Typeable, Data, Generic)
instance Show ChiSquared where
showsPrec i (ChiSquared n) = defaultShow1 "chiSquared" n i
instance Read ChiSquared where
readPrec = defaultReadPrecM1 "chiSquared" chiSquaredE
instance ToJSON ChiSquared
instance FromJSON ChiSquared where
parseJSON (Object v) = do
n <- v .: "chiSquaredNDF"
maybe (fail $ errMsg n) return $ chiSquaredE n
parseJSON _ = empty
instance Binary ChiSquared where
put (ChiSquared x) = put x
get = do n <- get
maybe (fail $ errMsg n) return $ chiSquaredE n
-- | Construct chi-squared distribution. Number of degrees of freedom
-- must be positive.
chiSquared :: Int -> ChiSquared
chiSquared n = maybe (error $ errMsg n) id $ chiSquaredE n
-- | Construct chi-squared distribution. Number of degrees of freedom
-- must be positive.
chiSquaredE :: Int -> Maybe ChiSquared
chiSquaredE n
| n <= 0 = Nothing
| otherwise = Just (ChiSquared n)
errMsg :: Int -> String
errMsg n = "Statistics.Distribution.ChiSquared.chiSquared: N.D.F. must be positive. Got " ++ show n
instance D.Distribution ChiSquared where
cumulative = cumulative
instance D.ContDistr ChiSquared where
density chi x
| x <= 0 = 0
| otherwise = exp $ log x * (ndf2 - 1) - x2 - logGamma ndf2 - log 2 * ndf2
where
ndf = fromIntegral $ chiSquaredNDF chi
ndf2 = ndf/2
x2 = x/2
logDensity chi x
| x <= 0 = m_neg_inf
| otherwise = log x * (ndf2 - 1) - x2 - logGamma ndf2 - log 2 * ndf2
where
ndf = fromIntegral $ chiSquaredNDF chi
ndf2 = ndf/2
x2 = x/2
quantile = quantile
instance D.Mean ChiSquared where
mean (ChiSquared ndf) = fromIntegral ndf
instance D.Variance ChiSquared where
variance (ChiSquared ndf) = fromIntegral (2*ndf)
instance D.MaybeMean ChiSquared where
maybeMean = Just . D.mean
instance D.MaybeVariance ChiSquared where
maybeStdDev = Just . D.stdDev
maybeVariance = Just . D.variance
instance D.Entropy ChiSquared where
entropy (ChiSquared ndf) =
let kHalf = 0.5 * fromIntegral ndf in
kHalf
+ log 2
+ logGamma kHalf
+ (1-kHalf) * digamma kHalf
instance D.MaybeEntropy ChiSquared where
maybeEntropy = Just . D.entropy
instance D.ContGen ChiSquared where
genContVar (ChiSquared n) = MWC.chiSquare n
cumulative :: ChiSquared -> Double -> Double
cumulative chi x
| x <= 0 = 0
| otherwise = incompleteGamma (ndf/2) (x/2)
where
ndf = fromIntegral $ chiSquaredNDF chi
quantile :: ChiSquared -> Double -> Double
quantile (ChiSquared ndf) p
| p == 0 = 0
| p == 1 = 1/0
| p > 0 && p < 1 = 2 * invIncompleteGamma (fromIntegral ndf / 2) p
| otherwise =
error $ "Statistics.Distribution.ChiSquared.quantile: p must be in [0,1] range. Got: "++show p
|