Calculate some interesting statistics about how hard it is to prove that dice aren't the normal distribution

distributions.hs

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]