{-# LANGUAGE DeriveDataTypeable #-}
-- |
-- Module    : Statistics.Distribution.LogNormal
-- Copyright : (c) 2009 Karamaan Group
--
-- The lognormal distribution. This is the distribution of a random 
-- variable whose logarithm is normally distributed.

module Statistics.Distribution.LogNormal
    (
      LogNormalDistribution
    -- * Constructors
    , fromParams
    , standard
    ) where

import Control.Exception (assert)
import Data.Number.Erf (erf)
import Data.Generics
import Statistics.Constants (m_sqrt_2, m_sqrt_2_pi)
import qualified Statistics.Distribution as D

-- | The lognormal distribution.
data LogNormalDistribution = ND {
      mean     :: {-# UNPACK #-} !Double
    , variance :: {-# UNPACK #-} !Double
    , ndPdfDenom :: {-# UNPACK #-} !Double
    , ndCdfDenom :: {-# UNPACK #-} !Double
    } deriving (Eq, Read, Show, Typeable, Data)

instance D.Distribution LogNormalDistribution where
    density    = density
    cumulative = cumulative
    quantile   = quantile

instance D.Variance LogNormalDistribution where
    variance = variance

instance D.Mean LogNormalDistribution where
    mean = mean

standard :: LogNormalDistribution
standard = ND {
             mean = 0.0
           , variance = 1.0
           , ndPdfDenom = m_sqrt_2_pi
           , ndCdfDenom = m_sqrt_2
           }

fromParams :: Double -> Double -> LogNormalDistribution
fromParams m v = assert (v > 0)
                 ND {
                   mean = m
                 , variance = v
                 , ndPdfDenom = m_sqrt_2_pi * sv
                 , ndCdfDenom = m_sqrt_2 * sv
                 }
    where sv = sqrt v

density :: LogNormalDistribution -> Double -> Double
density d x = exp (-xm * xm / (2 * variance d)) / (x * ndPdfDenom d)
    where xm = log x - mean d

cumulative :: LogNormalDistribution -> Double -> Double
cumulative d x = (1 + erf ((log x-mean d) / ndCdfDenom d)) / 2

-- | This is the quantile function for the LogNormalDistribution.
quantile :: LogNormalDistribution -> Double -> Double
quantile d p = exp $ quantile' d p

-- | This is the quantile function for NormalDistribution.
quantile' :: LogNormalDistribution -> Double -> Double
quantile' d p
  | p < 0 || p > 1 = inf/inf
  | p == 0         = -inf
  | p == 1         = inf
  | p == 0.5       = mean d
  | otherwise      = x * sqrt (variance d) + mean d
  where x          = D.findRoot standard p 0 (-100) 100
        inf        = 1/0