-
Notifications
You must be signed in to change notification settings - Fork 68
/
Normal.hs
131 lines (110 loc) · 4.16 KB
/
Normal.hs
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
{-# LANGUAGE BangPatterns, DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module : Statistics.Distribution.Normal
-- Copyright : (c) 2009 Bryan O'Sullivan
-- License : BSD3
--
-- Maintainer : bos@serpentine.com
-- Stability : experimental
-- Portability : portable
--
-- The normal distribution. This is a continuous probability
-- distribution that describes data that cluster around a mean.
module Statistics.Distribution.Normal
(
NormalDistribution
-- * Constructors
, normalDistr
, normalFromSample
, standard
) where
import Control.Applicative ((<$>), (<*>))
import Data.Binary (Binary)
import Data.Binary (put, get)
import Data.Data (Data, Typeable)
import GHC.Generics (Generic)
import Numeric.MathFunctions.Constants (m_sqrt_2, m_sqrt_2_pi)
import Numeric.SpecFunctions (erfc, invErfc)
import qualified Statistics.Distribution as D
import qualified Statistics.Sample as S
import qualified System.Random.MWC.Distributions as MWC
-- | The normal distribution.
data NormalDistribution = ND {
mean :: {-# UNPACK #-} !Double
, stdDev :: {-# UNPACK #-} !Double
, ndPdfDenom :: {-# UNPACK #-} !Double
, ndCdfDenom :: {-# UNPACK #-} !Double
} deriving (Eq, Read, Show, Typeable, Data, Generic)
instance Binary NormalDistribution where
put (ND w x y z) = put w >> put x >> put y >> put z
get = ND <$> get <*> get <*> get <*> get
instance D.Distribution NormalDistribution where
cumulative = cumulative
complCumulative = complCumulative
instance D.ContDistr NormalDistribution where
logDensity = logDensity
quantile = quantile
instance D.MaybeMean NormalDistribution where
maybeMean = Just . D.mean
instance D.Mean NormalDistribution where
mean = mean
instance D.MaybeVariance NormalDistribution where
maybeStdDev = Just . D.stdDev
maybeVariance = Just . D.variance
instance D.Variance NormalDistribution where
stdDev = stdDev
instance D.Entropy NormalDistribution where
entropy d = 0.5 * log (2 * pi * exp 1 * D.variance d)
instance D.MaybeEntropy NormalDistribution where
maybeEntropy = Just . D.entropy
instance D.ContGen NormalDistribution where
genContVar d = MWC.normal (mean d) (stdDev d)
{-# INLINE genContVar #-}
-- | Standard normal distribution with mean equal to 0 and variance equal to 1
standard :: NormalDistribution
standard = ND { mean = 0.0
, stdDev = 1.0
, ndPdfDenom = log m_sqrt_2_pi
, ndCdfDenom = m_sqrt_2
}
-- | Create normal distribution from parameters.
--
-- IMPORTANT: prior to 0.10 release second parameter was variance not
-- standard deviation.
normalDistr :: Double -- ^ Mean of distribution
-> Double -- ^ Standard deviation of distribution
-> NormalDistribution
normalDistr m sd
| sd > 0 = ND { mean = m
, stdDev = sd
, ndPdfDenom = log $ m_sqrt_2_pi * sd
, ndCdfDenom = m_sqrt_2 * sd
}
| otherwise =
error $ "Statistics.Distribution.Normal.normalDistr: standard deviation must be positive. Got " ++ show sd
-- | Create distribution using parameters estimated from
-- sample. Variance is estimated using maximum likelihood method
-- (biased estimation).
normalFromSample :: S.Sample -> NormalDistribution
normalFromSample xs
= normalDistr m (sqrt v)
where
(m,v) = S.meanVariance xs
logDensity :: NormalDistribution -> Double -> Double
logDensity d x = (-xm * xm / (2 * sd * sd)) - ndPdfDenom d
where xm = x - mean d
sd = stdDev d
cumulative :: NormalDistribution -> Double -> Double
cumulative d x = erfc ((mean d - x) / ndCdfDenom d) / 2
complCumulative :: NormalDistribution -> Double -> Double
complCumulative d x = erfc ((x - mean d) / ndCdfDenom d) / 2
quantile :: NormalDistribution -> Double -> Double
quantile d p
| p == 0 = -inf
| p == 1 = inf
| p == 0.5 = mean d
| p > 0 && p < 1 = x * ndCdfDenom d + mean d
| otherwise =
error $ "Statistics.Distribution.Normal.quantile: p must be in [0,1] range. Got: "++show p
where x = invErfc $ 2 * (1 - p)
inf = 1/0