gistfile1.hs

module Simpson where 

data Value = Value
    { minRange :: Double
    , maxRange :: Double
    , dof :: Int
    , expectedValue :: Double
    } deriving (Eq)

type Error = Double
type NumSegm = Double

simpsonRec :: Value -> NumSegm -> Error -> Double
simpsonRec val numSegm errorVal = if isValid val1 val2 then val1 else simpsonRec val (numSegm * 2) errorVal
                               where val1 = simpson val numSegm
                                     val2 = simpson val numSegm * 2
                                     isValid val1 val2 = (abs val1 val2) < errorVal

simpson :: Value -> NumSegm -> Double -> Double
simpson val numSegm errorVal = sum $ map (calculateTDistribution val w numSegm) list 
    where
    list = [0..numSegm]
    w =  (maxRange val) / numSegm

calculateTDistribution :: Value -> Double -> NumSegm -> Integer -> Double
calculateTDistribution values w numSegm x = (tDistribution((xi x) values) * multiplier(x numSegm)) * (w / 3)
                                            where xi x = w * x
multiplier i numSegm 
    | i == 0 || i == numSegm = 1
    | mod(i 2) == 0 = 2
    | otherwise = 4 

tDistribution xi values = gamma(dof + 1, 2) / exp(dof * pi, 0.5) * gamma(dof, 2)
                            * exp(1 + (square(xi) / dof), ((dof + 1) / 2) * (-1))

square = exp 2

data Fraction = Fraction Int Int
    deriving (Eq, Show)

gamma :: Integer -> Integer -> Double
gamma num den = if (mod num den) == 0 
                then factorial (floor(num / den)) - 1
                else gammaFrazione subtract((Fraction num den) 1)


factorial n = if n == 1 then 1 else n * factorial (n - 1)
oneHalf = Fraction 1 2
gammaFrazione :: Fraction -> Double
gammaFrazione f = if f == oneHalf 
                  then doubleValue oneHalf * sqrt(pi)
                  else doubleValue f * gammaFrazione $ subtractFrac f 1 

doubleValue (Fraction num den) = num / den
subtractFrac (Fraction num den) val = Fraction (num + (den * val)) den