{-# LANGUAGE CPP #-} {-# LANGUAGE ViewPatterns #-} import Control.Monad import Data.Word import Data.Proxy import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Unboxed.Mutable as MVU import Numeric.SpecFunctions (logChoose,incompleteGamma,log1p,logFactorial) import Test.Tasty import Test.Tasty.QuickCheck import Test.Tasty.Runners import Test.Tasty.Options import Test.Tasty.HUnit import Test.QuickCheck.Monadic import System.Random.MWC import System.Random.MWC.Distributions import System.Random.Stateful (StatefulGen) #if MIN_VERSION_random(1,3,0) import qualified System.Random.Stateful as Random (SeedGen(..)) #endif ---------------------------------------------------------------- -- ---------------------------------------------------------------- -- | Average number of events per bin for binned statistical tests newtype NPerBin = NPerBin Int instance IsOption NPerBin where defaultValue = NPerBin 100 parseValue = fmap NPerBin . safeRead optionName = pure "n-per-bin" optionHelp = pure "Average number of events per bin" -- | P-value for statistical test. newtype PValue = PValue Double instance IsOption PValue where defaultValue = PValue 1e-9 parseValue = fmap PValue . safeRead optionName = pure "pvalue" optionHelp = pure "P-value for statistical test" ---------------------------------------------------------------- -- ---------------------------------------------------------------- main :: IO () main = do -- Set up tasty let tasty_opts = [ Option (Proxy :: Proxy NPerBin) , Option (Proxy :: Proxy PValue) , Option (Proxy :: Proxy QuickCheckTests) , Option (Proxy :: Proxy QuickCheckReplay) , Option (Proxy :: Proxy QuickCheckShowReplay) , Option (Proxy :: Proxy QuickCheckMaxSize) , Option (Proxy :: Proxy QuickCheckMaxRatio) , Option (Proxy :: Proxy QuickCheckVerbose) , Option (Proxy :: Proxy QuickCheckMaxShrinks) ] ingredients = includingOptions tasty_opts : defaultIngredients opts <- parseOptions ingredients (testCase "Fake" (pure ())) let n_per_bin = lookupOption opts :: NPerBin p_val = lookupOption opts -- g0 <- createSystemRandom defaultMainWithIngredients ingredients $ testGroup "mwc" [ testProperty "save/restore" $ prop_SeedSaveRestore g0 #if MIN_VERSION_random(1,3,0) , testProperty "SeedGen" $ prop_SeedGen g0 #endif , testCase "user save/restore" $ saveRestoreUserSeed , testCase "empty seed data" $ emptySeed , testCase "output correct" $ do g <- create xs <- replicateM 513 (uniform g) assertEqual "[Word32]" xs golden , testCase "beta binomial mean" $ prop_betaBinomialMean , testProperty "binomial is binomial" $ prop_binomial_PMF n_per_bin p_val g0 , testProperty "Poisson is Poisson" $ prop_Poisson_PMF n_per_bin p_val g0 , testCase "poisson mean and variance of 1" $ prop_poissonMeanAndVar 1 , testCase "poisson mean and variance of 5" $ prop_poissonMeanAndVar 5 , testCase "poisson mean and variance of 40" $ prop_poissonMeanAndVar 40 , testCase "poisson mean and variance of 150" $ prop_poissonMeanAndVar 150 , testCase "poisson 0" $ do n <- poisson 0 g0 0 @=? n ] updateGenState :: GenIO -> IO () updateGenState g = replicateM_ 250 (uniform g :: IO Word32) prop_SeedSaveRestore :: GenIO -> Property prop_SeedSaveRestore g = monadicIO $ do run $ updateGenState g seed <- run $ save g seed' <- run $ save =<< restore seed return $ seed == seed' #if MIN_VERSION_random(1,3,0) prop_SeedGen :: GenIO -> Property prop_SeedGen g = monadicIO $ do run $ updateGenState g seed <- run $ save g return $ seed == (Random.fromSeed . Random.toSeed) seed #endif saveRestoreUserSeed :: IO () saveRestoreUserSeed = do let seed = toSeed $ U.replicate 258 0 seed' <- save =<< restore seed assertEqual "Seeds must be equal" seed' seed emptySeed :: IO () emptySeed = do let seed = toSeed U.empty seed' <- save =<< create assertEqual "Seeds must be equal" seed' seed -- First 513 values generated from seed made using create golden :: [Word32] golden = [ 2254043345, 562229898, 1034503294, 2470032534, 2831944869, 3042560015, 838672965, 715056843 , 3122641307, 2300516242, 4079538318, 3722020688, 98524204, 1450170923, 2669500465, 2890402829 , 114212910, 1914313000, 2389251496, 116282477, 1771812561, 1606473512, 1086420617, 3652430775 , 1165083752, 3599954795, 3006722175, 341614641, 3000394300, 1378097585, 1551512487, 81211762 , 604209599, 3949866361, 77745071, 3170410267, 752447516, 1213023833, 1624321744, 3251868348 , 1584957570, 2296897736, 3305840056, 1158966242, 2458014362, 1919777052, 3203159823, 3230279656 , 755741068, 3005087942, 2478156967, 410224731, 1196248614, 3302310440, 3295868805, 108051054 , 1010042411, 2725695484, 2201528637, 667561409, 79601486, 50029770, 566202616, 3217300833 , 2162817014, 925506837, 1527015413, 3079491438, 927252446, 118306579, 499811870, 2973454232 , 2979271640, 4078978924, 1864075883, 197741457, 296365782, 1784247291, 236572186, 464208268 , 1769568958, 827682258, 4247376295, 2959098022, 1183860331, 2475064236, 3952901213, 1953014945 , 393081236, 1616500498, 2201176136, 1663813362, 2167124739, 630903810, 113470040, 924745892 , 1081531735, 4039388931, 4118728223, 107819176, 2212875141, 1941653033, 3660517172, 192973521 , 3653156164, 1878601439, 3028195526, 2545631291, 3882334975, 456082861, 2775938704, 3813508885 , 1758481462, 3332769695, 3595846251, 3745876981, 152488869, 2555728588, 3058747945, 39382408 , 520595021, 2185388418, 3502636573, 2650173199, 1077668433, 3548643646, 71562049, 2726649517 , 494210825, 1208915815, 620990806, 2877290965, 3253243521, 804166732, 2481889113, 623399529 , 44880343, 183645859, 3283683418, 2214754452, 419328482, 4224066437, 1102669380, 1997964721 , 2437245376, 985749802, 858381069, 116806511, 1771295365, 97352549, 341972923, 2971905841 , 110707773, 950500868, 1237119233, 691764764, 896381812, 1528998276, 1269357470, 2567094423 , 52141189, 2722993417, 80628658, 3919817965, 3615946076, 899371181, 46940285, 4010779728 , 318101834, 30736609, 3577200709, 971882724, 1478800972, 3769640027, 3706909300, 3300631811 , 4057825972, 4285058790, 2329759553, 2967563409, 4080096760, 2762613004, 2518395275, 295718526 , 598435593, 2385852565, 2608425408, 604857293, 2246982455, 919156819, 1721573814, 2502545603 , 643962859, 587823425, 3508582012, 1777595823, 4119929334, 2833342174, 414044876, 2469473258 , 289159600, 3715175415, 966867024, 788102818, 3197534326, 3571396978, 3508903890, 570753009 , 4273926277, 3301521986, 1411959102, 2766249515, 4071012597, 959442028, 1962463990, 1098904190 , 714719899, 562204808, 1658783410, 1471669042, 2565780129, 1616648894, 4236521717, 1788863789 , 3068674883, 191936470, 253084644, 1915647866, 276372665, 2117183118, 3704675319, 218791054 , 3680045802, 406662689, 3844864229, 91140313, 3834015630, 25116147, 904830493, 3152559113 , 820358622, 1301896358, 296152699, 2202014455, 4256659428, 1175171414, 3287520873, 2028006499 , 327448717, 2095642873, 3798661296, 58567008, 3907537112, 3691259011, 1730142328, 2373011713 , 3387040741, 3189417655, 2949233059, 1238379614, 1813238023, 1064726446, 1339055235, 1744523609 , 279811576, 2934103599, 283542302, 994488448, 418691747, 1062780152, 102211875, 4071713296 , 1790834038, 1035092527, 2374272359, 3558280982, 1927663822, 3645417844, 3481790745, 3566282546 , 2000290859, 505518126, 363501589, 4075468679, 3247300709, 3705242654, 2731103609, 2836871038 , 589640144, 2546495106, 84767518, 1376911639, 2400770705, 527489676, 3804134352, 150084021 , 240070593, 3807594859, 3518576690, 659503830, 2239678479, 1273668921, 4271050554, 3090482972 , 401956859, 1772128561, 4438455, 1989666158, 2521484677, 3960178700, 4220196277, 1033999035 , 2214785840, 3428469341, 428564336, 2517446784, 3935757188, 3294001677, 1037971963, 3590324170 , 1220969729, 1719719817, 807688972, 77076422, 4251553858, 3963852375, 326128795, 3277818295 , 3671513069, 549617771, 1683950556, 3352913781, 409318429, 2456264774, 4036950639, 1162718475 , 83888874, 5578966, 172866494, 1542278848, 455546979, 1296511553, 4263636440, 2450589064 , 372411483, 211216338, 2632256495, 2393754408, 1336054289, 4087203071, 3159642437, 1933346856 , 2914152714, 3805541979, 2769740793, 1161287028, 2289749561, 4124509890, 2128452935, 210531695 , 4250709834, 390950534, 1421430300, 3030519715, 3228987297, 3086837053, 2866915453, 2335948692 , 1684378991, 2575634059, 4153427304, 2426048796, 4197556954, 2605152326, 2909410733, 2424889219 , 654577921, 811955499, 118126602, 504071559, 1278756230, 3896458168, 4105558075, 750276169 , 1120805572, 1762689330, 993728154, 1104363215, 774344996, 4077568952, 2183487324, 994724370 , 3323036885, 3880704963, 746305447, 961608310, 2030117337, 453935768, 800490463, 1034636 , 2323633564, 602565693, 806061242, 1899269713, 162686347, 467541008, 1529175313, 282891502 , 2529616339, 2930657178, 464272784, 2878535316, 807165854, 3209080518, 4080120278, 347748171 , 3972126063, 284174728, 2498328933, 1723872460, 143845955, 4223866687, 1761495357, 1544646770 , 4206103283, 3771574626, 642165282, 1119501013, 3514063332, 1443320304, 4056369796, 3602131475 , 1422908288, 804093687, 431176780, 40108717, 2998264213, 3705835674, 169805085, 454593842 , 2781536994, 2385225212, 4137367775, 2631435125, 2347082354, 629238010, 3283635219, 3815791831 , 1340400558, 4061846985, 3803921868, 3196119096, 718610843, 3694290834, 2169960411, 2407155570 , 2557480499, 16164105, 480957288, 2155919829, 2490067282, 2356287132, 511737296, 1602800634 , 1802275249, 3316832299, 50286484, 2106622541, 2352302834, 2538374315, 344766394, 2777260569 , 1215135803, 2229011963, 114632277, 1645499402, 1111617833, 3833259754, 928611385, 686744723 , 1898396834, 2461932251, 2665457318, 3797019621, 868313114, 2366635205, 481934875, 1170532970 , 642610859, 3150733309, 3508548582, 666714469, 711663449, 2436617656, 2681476315, 1637296693 , 2487349478, 4174144946, 2793869557, 559398604, 1898140528, 991962870, 864792875, 3861665129 , 4024051364, 3383200293, 773730975, 33517291, 2660126073, 689133464, 2248134097, 3874737781 , 3358012678] -- We can test two for the price of one betaBinomial :: StatefulGen g m => Double -> Double -> Int -> g -> m Int betaBinomial a b n g = do p <- beta a b g binomial n p g nSamples :: Int nSamples = 10000 alpha, delta :: Double alpha = 600.0 delta = 400.0 nTrials :: Int nTrials = 10 prop_betaBinomialMean :: IO () prop_betaBinomialMean = do g <- create ss <- replicateM nSamples $ betaBinomial alpha delta nTrials g let m = fromIntegral (sum ss) / fromIntegral nSamples let x1 = fromIntegral nTrials * alpha / (alpha + delta) assertBool ("Mean is " ++ show x1 ++ " but estimated as " ++ show m) (abs (m - x1) < 0.001) prop_poissonMeanAndVar :: Double -> IO () prop_poissonMeanAndVar lambda = do gen <- create ss <- replicateM nSamples $ poisson lambda gen let m = fromIntegral (sum ss) / fromIntegral nSamples :: Double let v = (fromIntegral (sum (map (^ 2) ss)) / fromIntegral nSamples) - (m ** 2) assertBool ("True mean and var: " ++ show lambda ++ " but estimated as " ++ show (m, v) ++ " respectively") ((abs (lambda - m) / lambda < 0.025) && (abs (lambda - v) / lambda < 0.025)) -- Test that `binomial` really samples from binomial distribution. -- -- If we have binomial random variate with number of trials N and -- sample it M times. Then number of events with K successes is -- described by multinomial distribution and we can test whether -- experimental distribution is described using likelihood ratio test prop_binomial_PMF :: NPerBin -> PValue -> GenIO -> Property prop_binomial_PMF (NPerBin n_per_bin) (PValue p_val) g = property $ do p <- choose (0, 1.0) -- Success probability n_trial <- choose (2, 100) -- Number of trials in binomial distribution -- Number of binomial samples to generate let n_samples = n_trial * n_per_bin n_samples' = fromIntegral n_samples -- Compute number of outcomes pure $ ioProperty $ do hist <- do buf <- MVU.new (n_trial + 1) replicateM_ n_samples $ MVU.modify buf (+(1::Int)) =<< binomial n_trial p g U.unsafeFreeze buf -- Here we compute twice log of likelihood ratio. Alternative -- hypothesis is some distribution which fits data perfectly -- -- Asymtotically it's ditributed as χ² with n_trial-1 degrees of -- freedom let likelihood _ 0 = 0 likelihood k (fromIntegral -> n_obs) = n_obs * (log (n_obs / n_samples') - logProbBinomial n_trial p k) let logL = 2 * U.sum (U.imap likelihood hist) let significance = 1 - cumulativeChi2 (n_trial - 1) logL pure $ counterexample ("p = " ++ show p) $ counterexample ("N = " ++ show n_trial) $ counterexample ("p-val = " ++ show significance) $ counterexample ("chi2 = " ++ show logL) $ significance > p_val -- Similar test for Poisson distribution. We build histogram by moving -- all values >λ+2σ into overflow bin. prop_Poisson_PMF :: NPerBin -> PValue -> GenIO -> Property prop_Poisson_PMF (NPerBin n_per_bin) (PValue p_val) g = property $ do lam <- choose (0, 100.0) -- Poisson rate let max_n = ceiling $ lam + 2 * sqrt lam -- Maximum histogram bin n_samples = max_n * n_per_bin -- Number of samples to generate n_samples' = fromIntegral n_samples pure $ ioProperty $ do hist <- do buf <- MVU.new (max_n + 1) replicateM_ n_samples $ do k <- poisson lam g MVU.modify buf (+(1::Int)) (min max_n k) U.unsafeFreeze buf let likelihood _ 0 = 0 likelihood k (fromIntegral -> n_obs) = n_obs * (log (n_obs / n_samples') - p) where p | k == max_n = logComlCdfPoisson lam k | otherwise = logProbPoisson lam k let logL = 2 * U.sum (U.imap likelihood hist) let significance = 1 - cumulativeChi2 max_n logL pure $ counterexample ("lambda = " ++ show lam) $ counterexample ("p-val = " ++ show significance) $ counterexample ("chi2 = " ++ show logL) $ significance > p_val ---------------------------------------------------------------- -- Statistical helpers ---------------------------------------------------------------- -- Logarithm of probability for binomial distribution logProbBinomial :: Int -> Double -> Int -> Double logProbBinomial n p k = logChoose n k + log p * k' + log1p (-p) * nk' where k' = fromIntegral k nk' = fromIntegral $ n - k logProbPoisson :: Double -> Int -> Double logProbPoisson lam i = log lam * fromIntegral i - logFactorial i - lam logComlCdfPoisson :: Double -> Int -> Double logComlCdfPoisson lam i = log $ incompleteGamma (fromIntegral (i :: Int)) lam cumulativeChi2 :: Int -> Double -> Double cumulativeChi2 (fromIntegral -> ndf) x | x <= 0 = 0 | otherwise = incompleteGamma (ndf/2) (x/2)