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
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
|
{-# OPTIONS_GHC -fno-warn-orphans #-}
{-# LANGUAGE ScopedTypeVariables #-}
-- Required for Param
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE OverlappingInstances #-}
{-# LANGUAGE ViewPatterns #-}
module Tests.Distribution (
distributionTests
) where
import Control.Applicative
import Control.Exception
import Data.List (find)
import Data.Typeable (Typeable)
import qualified Numeric.IEEE as IEEE
import Test.Framework (Test,testGroup)
import Test.Framework.Providers.QuickCheck2 (testProperty)
import Test.QuickCheck as QC
import Test.QuickCheck.Monadic as QC
import Text.Printf
import Statistics.Distribution
import Statistics.Distribution.Beta
import Statistics.Distribution.Binomial
import Statistics.Distribution.ChiSquared
import Statistics.Distribution.CauchyLorentz
import Statistics.Distribution.Exponential
import Statistics.Distribution.FDistribution
import Statistics.Distribution.Gamma
import Statistics.Distribution.Geometric
import Statistics.Distribution.Hypergeometric
import Statistics.Distribution.Normal
import Statistics.Distribution.Poisson
import Statistics.Distribution.StudentT
import Statistics.Distribution.Uniform
import Prelude hiding (catch)
import Tests.Helpers
-- | Tests for all distributions
distributionTests :: Test
distributionTests = testGroup "Tests for all distributions"
[ contDistrTests (T :: T BetaDistribution )
, contDistrTests (T :: T CauchyDistribution )
, contDistrTests (T :: T ChiSquared )
, contDistrTests (T :: T ExponentialDistribution )
, contDistrTests (T :: T GammaDistribution )
, contDistrTests (T :: T NormalDistribution )
, contDistrTests (T :: T UniformDistribution )
, contDistrTests (T :: T StudentT )
, contDistrTests (T :: T FDistribution )
, discreteDistrTests (T :: T BinomialDistribution )
, discreteDistrTests (T :: T GeometricDistribution )
, discreteDistrTests (T :: T HypergeometricDistribution )
, discreteDistrTests (T :: T PoissonDistribution )
, unitTests
]
----------------------------------------------------------------
-- Tests
----------------------------------------------------------------
-- Tests for continous distribution
contDistrTests :: (Param d, ContDistr d, QC.Arbitrary d, Typeable d, Show d) => T d -> Test
contDistrTests t = testGroup ("Tests for: " ++ typeName t) $
cdfTests t ++
[ testProperty "PDF sanity" $ pdfSanityCheck t
, testProperty "Quantile is CDF inverse" $ quantileIsInvCDF t
, testProperty "quantile fails p<0||p>1" $ quantileShouldFail t
]
-- Tests for discrete distribution
discreteDistrTests :: (Param d, DiscreteDistr d, QC.Arbitrary d, Typeable d, Show d) => T d -> Test
discreteDistrTests t = testGroup ("Tests for: " ++ typeName t) $
cdfTests t ++
[ testProperty "Prob. sanity" $ probSanityCheck t
, testProperty "CDF is sum of prob." $ discreteCDFcorrect t
]
-- Tests for distributions which have CDF
cdfTests :: (Param d, Distribution d, QC.Arbitrary d, Show d) => T d -> [Test]
cdfTests t =
[ testProperty "C.D.F. sanity" $ cdfSanityCheck t
, testProperty "CDF limit at +∞" $ cdfLimitAtPosInfinity t
, testProperty "CDF limit at -∞" $ cdfLimitAtNegInfinity t
, testProperty "CDF is nondecreasing" $ cdfIsNondecreasing t
, testProperty "1-CDF is correct" $ cdfComplementIsCorrect t
]
----------------------------------------------------------------
-- CDF is in [0,1] range
cdfSanityCheck :: (Distribution d) => T d -> d -> Double -> Bool
cdfSanityCheck _ d x = c >= 0 && c <= 1
where c = cumulative d x
-- CDF never decreases
cdfIsNondecreasing :: (Distribution d) => T d -> d -> Double -> Double -> Bool
cdfIsNondecreasing _ d = monotonicallyIncreasesIEEE $ cumulative d
-- CDF limit at +∞ is 1
cdfLimitAtPosInfinity :: (Param d, Distribution d) => T d -> d -> Property
cdfLimitAtPosInfinity _ d =
okForInfLimit d ==> printTestCase ("Last elements: " ++ show (drop 990 probs))
$ Just 1.0 == (find (>=1) probs)
where
probs = take 1000 $ map (cumulative d) $ iterate (*1.4) 1
-- CDF limit at -∞ is 0
cdfLimitAtNegInfinity :: (Param d, Distribution d) => T d -> d -> Property
cdfLimitAtNegInfinity _ d =
okForInfLimit d ==> printTestCase ("Last elements: " ++ show (drop 990 probs))
$ case find (< IEEE.epsilon) probs of
Nothing -> False
Just p -> p >= 0
where
probs = take 1000 $ map (cumulative d) $ iterate (*1.4) (-1)
-- CDF's complement is implemented correctly
cdfComplementIsCorrect :: (Distribution d) => T d -> d -> Double -> Bool
cdfComplementIsCorrect _ d x = (eq 1e-14) 1 (cumulative d x + complCumulative d x)
-- PDF is positive
pdfSanityCheck :: (ContDistr d) => T d -> d -> Double -> Bool
pdfSanityCheck _ d x = p >= 0
where p = density d x
-- Quantile is inverse of CDF
quantileIsInvCDF :: (Param d, ContDistr d) => T d -> d -> Double -> Property
quantileIsInvCDF _ d (snd . properFraction -> p) =
p > 0 && p < 1 ==> ( printTestCase (printf "Quantile = %g" q )
$ printTestCase (printf "Probability = %g" p )
$ printTestCase (printf "Probability' = %g" p')
$ printTestCase (printf "Error = %e" (abs $ p - p'))
$ abs (p - p') < invQuantilePrec d
)
where
q = quantile d p
p' = cumulative d q
-- Test that quantile fails if p<0 or p>1
quantileShouldFail :: (ContDistr d) => T d -> d -> Double -> Property
quantileShouldFail _ d p =
p < 0 || p > 1 ==> QC.monadicIO $ do r <- QC.run $ catch
(do { return $! quantile d p; return False })
(\(e :: SomeException) -> return True)
QC.assert r
-- Probability is in [0,1] range
probSanityCheck :: (DiscreteDistr d) => T d -> d -> Int -> Bool
probSanityCheck _ d x = p >= 0 && p <= 1
where p = probability d x
-- Check that discrete CDF is correct
discreteCDFcorrect :: (DiscreteDistr d) => T d -> d -> Int -> Int -> Property
discreteCDFcorrect _ d a b
= printTestCase (printf "CDF = %g" p1)
$ printTestCase (printf "Sum = %g" p2)
$ printTestCase (printf "Δ = %g" (abs (p1 - p2)))
$ abs (p1 - p2) < 3e-10
-- Avoid too large differeneces. Otherwise there is to much to sum
--
-- Absolute difference is used guard againist precision loss when
-- close values of CDF are subtracted
where
n = min a b
m = n + (abs (a - b) `mod` 100)
p1 = cumulative d (fromIntegral m + 0.5) - cumulative d (fromIntegral n - 0.5)
p2 = sum $ map (probability d) [n .. m]
----------------------------------------------------------------
-- Arbitrary instances for ditributions
----------------------------------------------------------------
instance QC.Arbitrary BinomialDistribution where
arbitrary = binomial <$> QC.choose (1,100) <*> QC.choose (0,1)
instance QC.Arbitrary ExponentialDistribution where
arbitrary = exponential <$> QC.choose (0,100)
instance QC.Arbitrary GammaDistribution where
arbitrary = gammaDistr <$> QC.choose (0.1,10) <*> QC.choose (0.1,10)
instance QC.Arbitrary BetaDistribution where
arbitrary = betaDistr <$> QC.choose (1e-3,10) <*> QC.choose (1e-3,10)
instance QC.Arbitrary GeometricDistribution where
arbitrary = geometric <$> QC.choose (0,1)
instance QC.Arbitrary HypergeometricDistribution where
arbitrary = do l <- QC.choose (1,20)
m <- QC.choose (0,l)
k <- QC.choose (1,l)
return $ hypergeometric m l k
instance QC.Arbitrary NormalDistribution where
arbitrary = normalDistr <$> QC.choose (-100,100) <*> QC.choose (1e-3, 1e3)
instance QC.Arbitrary PoissonDistribution where
arbitrary = poisson <$> QC.choose (0,1)
instance QC.Arbitrary ChiSquared where
arbitrary = chiSquared <$> QC.choose (1,100)
instance QC.Arbitrary UniformDistribution where
arbitrary = do a <- QC.arbitrary
b <- QC.arbitrary `suchThat` (/= a)
return $ uniformDistr a b
instance QC.Arbitrary CauchyDistribution where
arbitrary = cauchyDistribution
<$> arbitrary
<*> ((abs <$> arbitrary) `suchThat` (> 0))
instance QC.Arbitrary StudentT where
arbitrary = studentT <$> ((abs <$> arbitrary) `suchThat` (>0))
instance QC.Arbitrary FDistribution where
arbitrary = fDistribution
<$> ((abs <$> arbitrary) `suchThat` (>0))
<*> ((abs <$> arbitrary) `suchThat` (>0))
-- Parameters for distribution testing. Some distribution require
-- relaxing parameters a bit
class Param a where
-- Precision for quantileIsInvCDF
invQuantilePrec :: a -> Double
invQuantilePrec _ = 1e-14
-- Distribution is OK for testing limits
okForInfLimit :: a -> Bool
okForInfLimit _ = True
instance Param a
instance Param StudentT where
invQuantilePrec _ = 1e-13
okForInfLimit d = studentTndf d > 0.75
instance Param FDistribution where
invQuantilePrec _ = 1e-12
----------------------------------------------------------------
-- Unit tests
----------------------------------------------------------------
unitTests :: Test
unitTests = testGroup "Unit tests"
[ testAssertion "density (gammaDistr 150 1/150) 1 == 4.883311" $
4.883311418525483 =~ (density (gammaDistr 150 (1/150)) 1)
-- Student-T
, testStudentPDF 0.3 1.34 0.0648215 -- PDF
, testStudentPDF 1 0.42 0.27058
, testStudentPDF 4.4 0.33 0.352994
, testStudentCDF 0.3 3.34 0.757146 -- CDF
, testStudentCDF 1 0.42 0.626569
, testStudentCDF 4.4 0.33 0.621739
-- F-distribution
, testFdistrPDF 1 3 3 (1/(6 * pi)) -- PDF
, testFdistrPDF 2 2 1.2 0.206612
, testFdistrPDF 10 12 8 0.000385613179281892790166
, testFdistrCDF 1 3 3 0.81830988618379067153 -- CDF
, testFdistrCDF 2 2 1.2 0.545455
, testFdistrCDF 10 12 8 0.99935509863451408041
]
where
-- Student-T
testStudentPDF ndf x exact
= testAssertion (printf "density (studentT %f) %f ≈ %f" ndf x exact)
$ eq 1e-5 exact (density (studentT ndf) x)
testStudentCDF ndf x exact
= testAssertion (printf "cumulative (studentT %f) %f ≈ %f" ndf x exact)
$ eq 1e-5 exact (cumulative (studentT ndf) x)
-- F-distribution
testFdistrPDF n m x exact
= testAssertion (printf "density (fDistribution %i %i) %f ≈ %f [got %f]" n m x exact d)
$ eq 1e-5 exact d
where d = density (fDistribution n m) x
testFdistrCDF n m x exact
= testAssertion (printf "cumulative (fDistribution %i %i) %f ≈ %f [got %f]" n m x exact d)
$ eq 1e-5 exact d
where d = cumulative (fDistribution n m) x
|
