Created
April 10, 2015 18:18
-
-
Save Cedev/a0779f0666bdf02ddf62 to your computer and use it in GitHub Desktop.
Calculate some interesting statistics about how hard it is to prove that dice aren't the normal distribution
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import Statistics.Distribution | |
import Statistics.Distribution.Normal | |
data D n = D {n :: n, s :: n} | |
instance Integral n => Distribution (D n) where | |
cumulative d x = cumulativeD d (floor x) | |
instance Integral n => DiscreteDistr (D n) where | |
probability d x = probabilityD d (fromIntegral x) | |
instance Integral n => MaybeMean (D n) where | |
maybeMean = Just . meanD | |
instance Integral n => Mean (D n) where | |
mean = meanD | |
meanD :: (Integral n, Fractional x) => D n -> x | |
meanD (D n s) = fromIntegral n * (fromIntegral s + 1)/2 | |
varD :: (Integral n, Fractional x) => D n -> x | |
varD (D 1 s) = (fromIntegral s ^ 2 - 1)/12 | |
varD (D n s) = fromIntegral n * varD (D 1 s) | |
rangeD :: (Integral n) => D n -> [n] | |
rangeD (D n s) = [n .. n * s] | |
-- from http://mathworld.wolfram.com/Dice.html | |
coefD :: (Integral n) => D n -> n -> n | |
coefD (D n s) p = sum . map f $ [0 .. (p - n) `div` s] | |
where | |
f k = ((-1)^k) * (n `choose` k) * ((p - (s * k) - 1) `choose` (n - 1)) | |
coefsD d = table (coefD d) (rangeD d) | |
probabilityD :: (Integral n, Fractional x) => D n -> n -> x | |
probabilityD d x = fromIntegral (coefD d x)/fromIntegral (s d^n d) | |
densitiesD :: (Integral n, Fractional x) => D n -> [(n, x)] | |
densitiesD d = table (probabilityD d) (rangeD d) | |
cumulativeD :: (Integral n, Fractional x) => D n -> n -> x | |
cumulativeD d x = sum . map snd . takeWhile ((<= x) . fst) . densitiesD $ d | |
distrD :: (Integral n) => D n -> NormalDistribution | |
distrD d = normalDistr (meanD d) (sqrt (varD d)) | |
normalPsD d = table (cumulative nd . fromIntegral) (rangeD d) | |
where | |
nd = distrD d | |
cumsD d = table (cumulative d . fromIntegral) (rangeD d) | |
choose :: (Integral n) => n -> n -> n | |
choose n 0 = 1 | |
choose 0 k = 0 | |
choose n k = choose (n-1) (k-1) * n `div` k | |
table :: (a -> b) -> [a] -> [(a, b)] | |
table f = map (\a -> (a, f a)) | |
-- from http://stats.stackexchange.com/q/145583/35181 | |
nChooseBernoulli :: Floating a => a -> a -> a -> a | |
nChooseBernoulli z q1 q2 = (z/(phi q2 - phi q1))^2 | |
where | |
phi = asin . sqrt | |
{- | |
nShowNotNormalD :: (Integral n, Real z) => z -> D n -> [(Double, Double)] | |
nShowNotNormalD z d@(D n s) = table f . map ((+(5/10)) . fromIntegral) $ (n - 1 : rangeD d) | |
where | |
f :: Double -> Double | |
f x = nChooseBernoulli (realToFrac z) (cumulative d x) (cumulative nd x) | |
nd = distrD d | |
-} | |
printTable :: (Integral n) => D n -> IO () | |
printTable d@(D n s) = do | |
print $ "D " ++ show (toInteger n) ++ " " ++ show (toInteger s) | |
let points = [o + fromIntegral x | x <- (rangeD d), o <- [0]] | |
nd = distrD d | |
printLine p = do | |
putStr . show . round $ p | |
putStr "\t" | |
putStr . show . cumulative d $ p | |
putStr "\t" | |
putStr . show . cumulative nd $ p | |
putStr "\t" | |
putStr . show . ceiling $ nChooseBernoulli 1 (cumulative d p) (cumulative nd p) | |
--putStr "\t" | |
--putStr . show . ceiling $ nChooseBernoulli 1 (1 - cumulative d p) (1 - cumulative nd p) | |
putStr "\t" | |
putStr . show . ceiling $ nChooseBernoulli 2 (cumulative d p) (cumulative nd p) | |
putStrLn "" | |
mapM_ printLine points | |
main = do | |
let cups = [D 1 6, D 2 6, D 3 6, D 4 6, D 6 6, D 4 3, D 8 3, D 1 20] | |
let stats = [ | |
putStrLn . const "", | |
print . rangeD, | |
print . meanD, | |
print . varD, | |
print . coefsD, | |
print . cumsD, | |
print . distrD, | |
print . normalPsD, | |
-- print . nShowNotNormalD 1, | |
printTable | |
] | |
sequence_ [s c | c <- cups, s <- stats] |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment