    {-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
    -- |
    -- Module    : Statistics.Distribution.StudentT
    -- Copyright : (c) 2011 Aleksey Khudyakov
    -- License   : BSD3
    --
    -- Maintainer  : bos@serpentine.com
    -- Stability   : experimental
    -- Portability : portable
    --
    -- Student-T distribution
    module Statistics.Distribution.StudentT (
        StudentT
      , studentT
      , studentTndf
      , studentTUnstandardized
      ) where

    import Data.Aeson (FromJSON, ToJSON)
    import Data.Binary (Binary)
    import Data.Data (Data, Typeable)
    import GHC.Generics (Generic)
    import qualified Statistics.Distribution as D
    import Statistics.Distribution.Transform (LinearTransform (..))
    import Numeric.SpecFunctions (
      logBeta, incompleteBeta, invIncompleteBeta, digamma)
    import Data.Binary (put, get)

    -- | Student-T distribution
    newtype StudentT = StudentT { studentTndf :: Double }
                       deriving (Eq, Show, Read, Typeable, Data, Generic)

    instance FromJSON StudentT
    instance ToJSON StudentT

    instance Binary StudentT where
        put = put . studentTndf
        get = fmap StudentT get

    -- | Create Student-T distribution. Number of parameters must be positive.
    studentT :: Double -> StudentT
    studentT ndf
      | ndf > 0   = StudentT ndf
      | otherwise = modErr "studentT" "non-positive number of degrees of freedom"

    instance D.Distribution StudentT where
      cumulative = cumulative

    instance D.ContDistr StudentT where
      density    d@(StudentT ndf) x = exp (logDensityUnscaled d x) / sqrt ndf
      logDensity d@(StudentT ndf) x = logDensityUnscaled d x - log (sqrt ndf)
      quantile = quantile

    cumulative :: StudentT -> Double -> Double
    cumulative (StudentT ndf) x
      | x > 0     = 1 - 0.5 * ibeta
      | otherwise = 0.5 * ibeta
      where
        ibeta = incompleteBeta (0.5 * ndf) 0.5 (ndf / (ndf + x*x))

    logDensityUnscaled :: StudentT -> Double -> Double
    logDensityUnscaled (StudentT ndf) x =
        log (ndf / (ndf + x*x)) * (0.5 * (1 + ndf)) - logBeta 0.5 (0.5 * ndf)

    quantile :: StudentT -> Double -> Double
    quantile (StudentT ndf) p
      | p >= 0 && p <= 1 =
        let x = invIncompleteBeta (0.5 * ndf) 0.5 (2 * min p (1 - p))
        in case sqrt $ ndf * (1 - x) / x of
             r | p < 0.5   -> -r
               | otherwise -> r
      | otherwise = modErr "quantile" $ "p must be in [0,1] range. Got: "++show p

    instance D.MaybeMean StudentT where
      maybeMean (StudentT ndf) | ndf > 1   = Just 0
                               | otherwise = Nothing

    instance D.MaybeVariance StudentT where
      maybeVariance (StudentT ndf) | ndf > 2   = Just $! ndf / (ndf - 2)
                                   | otherwise = Nothing

    instance D.Entropy StudentT where
      entropy (StudentT ndf) =
        0.5 * (ndf+1) * (digamma ((1+ndf)/2) - digamma(ndf/2))
        + log (sqrt ndf)
        + logBeta (ndf/2) 0.5

    instance D.MaybeEntropy StudentT where
      maybeEntropy = Just . D.entropy

    instance D.ContGen StudentT where
      genContVar = D.genContinous

    -- | Create an unstandardized Student-t distribution.
    studentTUnstandardized :: Double -- ^ Number of degrees of freedom
                           -> Double -- ^ Central value (0 for standard Student T distribution)
                           -> Double -- ^ Scale parameter
                           -> LinearTransform StudentT
    studentTUnstandardized ndf mu sigma
      | sigma > 0 = LinearTransform mu sigma $ studentT ndf
      | otherwise = modErr "studentTUnstandardized" $ "sigma must be > 0. Got: " ++ show sigma

    modErr :: String -> String -> a
    modErr fun msg = error $ "Statistics.Distribution.StudentT." ++ fun ++ ": " ++ msg