ベルヌイ試行・二項分布

BernoulliTrial.hs

module Bernoulli (
    f,
    p,
    c,
    bp,
    be,
    bv
) where

{--
    calcurate the fractorical of non-negative integer
    @param n
    @return n!
--}
f :: (Integral a) => (Real b) => a -> b 
f n
    | 0 > n = error "n must be >=0."
f 0 = 1
f n = n' * f (n-1)
    where n' = fromIntegral n

{--
    calcurate the permutation.
    @param n
    @param i <= n
    @return nPi
--}
p :: (Integral a) => (Real b) => a -> a -> b 
n`p`i
    | n < i     = error "n must be >=i."
_`p`0 = 1
n`p`i = n' * ((n-1)`p`(i-1))
    where n' = fromIntegral n

{--
    calcurate the combination. 
    @param n
    @param i <= n
    @return nCi
--}
c :: (Integral a) => (Real b) => a -> a -> b
n`c`i = fromIntegral $ (n`p`i) `div` (f i)

{--
    calculate the probability distribution for Bernulli Trial
    @param n number of experiments
    @param p probability of success experiment
    @param i number of success experiment
    @param the probability of having i success experiments after n experiments
--}
bp :: (Integral a) => (RealFloat b) =>
    a -> b -> a -> b
bp n p i
    | n < i            = error "n must be >i."
    | p < 0 || p > 1   = error "p must be 0< and 1>."
bp n p i = (n`c`i) * (p**i') * (q**(n'-i'))
    where q = 1 - p; n' = fromIntegral n; i' = fromIntegral i;

{--
    calcurate the expected value for Bernulli Trial
    @param n number of experiments
    @param p probability of success experiment
    @param the expected value of success experiments after n experiments
--}
be :: (Integral a) => (RealFloat b) => 
    a -> b -> b 
be n p = sum [(fromIntegral i) * bp n p i | i<-[0..n]]

{--
    calcurate the valiance for Bernulli Trial
    @param n number of experiments
    @param p probability of success experiment
    @param the valiance of success experiments after n experiments
--}
bv :: (Integral a) => (RealFloat b) =>
    a -> b -> b
bv n p = (sum [(fromIntegral i)**2 * bp n p i | i<-[0..n]]) - ((be n p)**2)

LaplacianLow.hs

import Bernoulli
import System.Environment

main = do
    args <- getArgs
    let
        n = read (args !! 0) :: Int
        l = read (args !! 1) :: Double
        p = l / (fromIntegral n)
    sequence [print $ Bernoulli.bp n p i | i <- [0..n]]