gev-lib (empty) → 0.1.0.0
raw patch · 8 files changed
+534/−0 lines, 8 filesdep +HUnitdep +basedep +gev-lib
Dependencies added: HUnit, base, gev-lib, random
Files
- CHANGELOG.md +5/−0
- gev-lib.cabal +45/−0
- src/Gev.hs +50/−0
- src/Gev/Frechet.hs +77/−0
- src/Gev/GevDist.hs +98/−0
- src/Gev/Gumbel.hs +68/−0
- src/Gev/Weibull.hs +72/−0
- src/Tests/tests.hs +119/−0
+ CHANGELOG.md view
@@ -0,0 +1,5 @@+# Revision history for gev-test++## 0.1.0.0 -- YYYY-mm-dd++* First version. Released on an unsuspecting world.
+ gev-lib.cabal view
@@ -0,0 +1,45 @@+cabal-version: 2.4+name: gev-lib+version: 0.1.0.0++-- A short (one-line) description of the package.+synopsis: The family of Extreme Value Distributions.++-- A longer description of the package.+description: Basic Distributional quantities (CDF, PDF, Quantile) for the Gumbel, Féchet, Weibull and GEV Distributions.++-- A URL where users can report bugs.+-- bug-reports:++-- The license under which the package is released.+license: ISC++-- The package author(s).+author: Gabriel Haeck++-- An email address to which users can send suggestions, bug reports, and patches.+maintainer: haeckgabriel@gmail.com++-- A copyright notice.+-- copyright:+category: Statistics+extra-source-files: CHANGELOG.md++library gev-dist+ exposed-modules: Gev, Gev.Gumbel, Gev.Frechet, Gev.Weibull, Gev.GevDist++ -- Modules included in this library but not exported.+ -- other-modules:++ -- LANGUAGE extensions used by modules in this package.+ -- other-extensions:+ build-depends: base ^>=4.14.3.0, random+ hs-source-dirs: src+ default-language: Haskell2010++test-suite tests+ type: exitcode-stdio-1.0+ hs-source-dirs: .+ main-is: src/Tests/tests.hs+ build-depends: base ^>=4.14.3.0, HUnit ^>=1.6, gev-dist+ default-language: Haskell2010
+ src/Gev.hs view
@@ -0,0 +1,50 @@+module Gev ( Distribution(..) ) where++import System.Random.Stateful (StatefulGen, uniformDouble01M)++-- | Distribution Class for the GEV family of distributions. That is, each of the +-- distributions considered will have a CDF, PDF and Quantile function.+class Distribution d where+ -- | Cumulative Distribution Function (CDF) of a given distribution.+ -- i.e. $\mathbb{P}(X \leq x)$ for $x \in \Omega(X)$ (i.e. x is in the support of X)+ --+ -- > cdf d +∞ = 1+ -- > cdf d -∞ = 0+ cdf :: d -> Double -> Double+ cdf d x = 1 - complCdf d x++ -- | Complement of the CDF, i.e. $\mathbb{P}(X \geq x)$.+ complCdf :: d -> Double -> Double+ complCdf d x = 1 - cdf d x++ -- | Probability Density Function (pdf) of a distribution.+ -- i.e. $\mathbb{P}(X = x)$ for $x \in \Omega(X)$ (i.e. x is in the support of X)+ pdf :: d -> Double -> Double+ pdf d = exp . logPdf d+ + -- | Log density of a given distribution+ -- i.e. density for $Y = \log X$+ logPdf :: d -> Double -> Double+ logPdf d = log . pdf d++ -- | Quantile function (a.k.a inverse CDF) of a distribution.+ -- i.e. $F^{-1}(x)$ for $x \in [0, 1]$.+ quantile :: d -> Double -> Double+ quantile d x = complQuantile d (1- x)++ -- | Quantile complement, i.e. Quantile for level $1 - \alpha$.+ complQuantile :: d -> Double -> Double+ complQuantile d x = quantile d (1 - x)++ -- | generate random value of the Distribution.+ randGen :: StatefulGen g m => d -> g -> m Double+ randGen d gen = do+ x <- uniformDouble01M gen+ return $! quantile d x+ {-# MINIMAL (cdf | complCdf), (pdf | logPdf), (quantile | complQuantile) #-}++-- see https://hackage.haskell.org/package/statistics-0.16.1.2/docs/src/Statistics.Distribution.html#genContinuous+--genContinuous :: (ContDistr d, StatefulGen g m) => d -> g -> m Double+--genContinuous d gen = do+-- x <- uniformDouble01M gen+-- return $! quantile d x
+ src/Gev/Frechet.hs view
@@ -0,0 +1,77 @@+module Gev.Frechet+ (+ FrechetDistribution+ -- * constructors+ , frechetDist+ , frechetDistMaybe+ -- * accessors+ , location+ , scale+ , shape+ ) where++import qualified Gev++data FrechetDistribution = Frechet {+ location :: {-# UNPACK #-} !Double+ , scale :: {-# UNPACK #-} !Double+ , shape :: {-# UNPACK #-} !Double+ } deriving (Eq)++instance Show FrechetDistribution where+ show (Frechet loc sc sh) = show "Frechet Distribution; loc: " ++ show loc ++ ", scale: " ++ show sc ++ " and shape: " ++ show sh++-- error message when initiating Frechet with scale parameter less than 0.+frechetErrMsg :: Double -> Double -> Double -> String+frechetErrMsg loc scale sh = "Gev.Frechet: " + ++ "loc = " ++ show loc+ ++ " scale = " ++ show scale+ ++ " schape = " ++ show sh+ ++ ", but both the scale and shape parameters must be positive!"++-- | create Frechet Dist, where scale parameter must be greater than 0.+frechetDistMaybe :: Double -> Double -> Double -> Maybe FrechetDistribution+frechetDistMaybe loc sc sh+ | sc > 0 && sh > 0 = Just $ Frechet loc sc sh+ | otherwise = Nothing++-- | create Frechet Dist, where scale parameter must be greater than 0.+frechetDist :: Double -> Double -> Double -> FrechetDistribution+frechetDist loc sc sh = maybe (error $ frechetErrMsg loc sc sh) id $ frechetDistMaybe loc sc sh++-- | The CDF of the Frechet distribution+cdfFrechet :: FrechetDistribution -> Double -> Double+cdfFrechet (Frechet loc sc sh) x+ | x <= 0 = 0+ | otherwise = exp $ - (y ** pow)+ where + y = (x - loc) / sc+ pow = - sh++-- | The PDF of the Frechet distribution+pdfFrechet :: FrechetDistribution -> Double -> Double+pdfFrechet (Frechet loc sc sh) x =+ let y = (x - loc) / sc+ const = sh / sc+ expterm = exp $ - (y ** (- sh))+ middle = y ** (-1 - sh)+ in const * middle * expterm++-- Quantile function of the Frechet Distribution+quantileFrechet :: FrechetDistribution -> Double -> Double+quantileFrechet (Frechet loc sc sh) x+ | x == 0 = -inf+ | x == 1 = inf+ | x > 0 && x < 1 = loc + sc * (logexp ** pow)+ | otherwise =+ error $ "Gev.FrechetDistribution.quantile: The given value must be between 0 and 1, got: " ++ show x+ where + inf = 1 / 0+ logexp = - log x+ pow = - 1 / sh++-- | Gev.Distribution instance implementation for the Frechet Distribution+instance Gev.Distribution FrechetDistribution where+ cdf = cdfFrechet+ pdf = pdfFrechet+ quantile = quantileFrechet
+ src/Gev/GevDist.hs view
@@ -0,0 +1,98 @@+module Gev.GevDist + (+ GevDistribution+ -- * constructors+ , gevDist+ , gevDistMaybe+ -- * accessors+ , location+ , scale+ , shape+ ) where++import qualified Gev++data GevDistribution = GEV {+ location :: {-# UNPACK #-} !Double+ , scale :: {-# UNPACK #-} !Double+ , shape :: {-# UNPACK #-} !Double+ } deriving (Eq)++instance Show GevDistribution where+ show (GEV loc sc sh) = show "GEV Distribution; loc: " ++ show loc ++ ", scale: " ++ show sc ++ " and shape: " ++ show sh++-- error message when initiating GEV Distribution with scale parameter less than 0.+gevErrMsg :: Double -> Double -> Double -> String+gevErrMsg loc scale sh = "Gev.GevDist: " + ++ "loc = " ++ show loc+ ++ " scale = " ++ show scale+ ++ " schape = " ++ show sh+ ++ ", but the scale parameter must be positive!"++gevDistMaybe :: Double -> Double -> Double -> Maybe GevDistribution+gevDistMaybe loc sc sh+ | sc > 0 = Just $ GEV loc sc sh+ | otherwise = Nothing++gevDist :: Double -> Double -> Double -> GevDistribution+gevDist loc sc sh = maybe (error $ gevErrMsg loc sc sh) id $ gevDistMaybe loc sh sc++-- | helper function: checks that the given x value falls in the valid part of the GEV domain.+-- -- actually this i isn't needed, just included as a guard pattern it'll be fine..+gevDomainCheck :: GevDistribution -> Double -> Bool+gevDomainCheck (GEV loc sc sh) x =+ if 1 + sh * (x - loc / sc) > 0 then True else False++-- | t(x) function depends on if the shape parameter (\zeta) is 0 or not.+-- t(x) = \exp \left(x) = \left( 1 + \zeta \left( \frac{x - \mu}{ \sigma} \right) \right)^{- \frac{1}{\zeta}}$$ if $\zeta \neq 0$,+-- or $t(x) = \exp \left \{ - \frac{x - \mu}{ \sigma } \right \}$ if $\zeta = 0$+gevArg :: Double -> Double -> Double -> Double -> Double+gevArg loc sc sh x =+ if sh == 0 then one else two+ where+ y = (x - loc) / sc+ one = exp $ - y+ two = (1 + sh * y) ** (- 1 / sh)+ +-- | CDF of the GEV Distribution.+cdfGEV :: GevDistribution -> Double -> Double+cdfGEV (GEV loc sc sh) x+ | x <= 0 = 0+ | 1 + sh * (x - loc / sc) > 0 = exp $ - tVal+ | otherwise =+ error $ "Gev.GevDist.cdf: The given x value is not in the support of the Distribution: " ++ show x+ where+ tVal = gevArg loc sc sh x+ ++-- | PDF of the GEV Distribution.+pdfGEV :: GevDistribution -> Double -> Double+pdfGEV (GEV loc sc sh) x+ | 1 + sh * (x - loc / sc) > 0 = const * middle * exp ( - tVal)+ | otherwise =+ error $ "Gev.GevDist.pdf: The given x value is not in the support of the Distribution: " ++ show x+ where+ tVal = gevArg loc sc sh x+ const = 1 / sc+ middle = tVal ** (sh + 1)++-- | Quantile function of the GEV Distribution+-- Quantile function of the Frechet Distribution+quantileGEV :: GevDistribution -> Double -> Double+quantileGEV (GEV loc sc sh) x+ | x == 0 = -inf+ | x == 1 = inf+ | x > 0 && x < 1 = if sh == 0 then one else two+ | otherwise =+ error $ "Gev.GEVDist.quantile: The given value must be between 0 and 1, got: " ++ show x+ where + inf = 1 / 0+ logx = - log x+ one = - sc * log logx + loc+ const = sc / sh+ two = const * (logx ** (- sh)) - const + loc++instance Gev.Distribution GevDistribution where+ cdf = cdfGEV+ pdf = pdfGEV+ quantile = quantileGEV
+ src/Gev/Gumbel.hs view
@@ -0,0 +1,68 @@+module Gev.Gumbel+ (+ GumbelDistribution+ -- * Constructors+ , gumbelDist+ , gumbelDistMaybe+ -- * Accessors+ , location+ , scale+ ) where++import qualified Gev+-- import Numeric.SpecFunctions (log1p,expm1) POSSIBLY!!!++data GumbelDistribution = Gumbel {+ location :: {-# UNPACK #-} !Double+ , scale :: {-# UNPACK #-} !Double+ } deriving (Eq)++instance Show GumbelDistribution where+ show (Gumbel loc sc) = show "Gumbel Distribution; loc: " ++ show loc ++ " and scale: " ++ show sc++-- error message when initiating Gumbel with scale parameter less than 0.+gumbelErrMsg :: Double -> Double -> String+gumbelErrMsg loc scale = "Gev.Gumbel: " + ++ "loc = " ++ show loc+ ++ " scale = " ++ show scale+ ++ ", but the scale parameter must be positive!"++-- | create Gumbel Dist, where scale parameter must be greater than 0.+gumbelDistMaybe :: Double -> Double -> Maybe GumbelDistribution+gumbelDistMaybe loc sc+ | sc > 0 = Just $ Gumbel loc sc+ | otherwise = Nothing++-- | create Gumbel Dist, where scale parameter must be greater than 0.+gumbelDist :: Double -> Double -> GumbelDistribution+gumbelDist loc sc = maybe (error $ gumbelErrMsg loc sc) id $ gumbelDistMaybe loc sc++-- | The CDF of the Gumbel distribution+cdfGumbel :: GumbelDistribution -> Double -> Double+cdfGumbel (Gumbel loc sc) x+ | x <= 0 = 0+ | otherwise = exp $ - (exp y)+ where + y = - (x - loc) / sc++-- | The PDf of the Gumbel distribution+pdfGumbel :: GumbelDistribution -> Double -> Double+pdfGumbel (Gumbel loc sc) x =+ let exp_y = exp $ - (x- loc) / sc+ const = 1 / sc+ in const * exp_y * (exp $ - exp_y)++quantileGumbel :: GumbelDistribution -> Double -> Double+quantileGumbel (Gumbel loc sc) x+ | x == 0 = -inf+ | x == 1 = inf+ | x > 0 && x < 1 = loc - sc * (log $ - log x)+ | otherwise =+ error $ "Gev.GumbelDistribution.quantile: The given value must be between 0 and 1, got: " ++ show x+ where inf = 1 / 0++-- | Gev.Distribution instance implementation for the Gumbel Distribution+instance Gev.Distribution GumbelDistribution where+ cdf = cdfGumbel+ pdf = pdfGumbel+ quantile = quantileGumbel
+ src/Gev/Weibull.hs view
@@ -0,0 +1,72 @@+module Gev.Weibull+ (+ WeibullDistribution+ -- * constructors+ , weibullDist+ , weibullDistMaybe+ -- * accessors+ , location+ , scale+ , shape+ ) where++import qualified Gev++data WeibullDistribution = Weibull {+ location :: {-# UNPACK #-} !Double+ , scale :: {-# UNPACK #-} !Double+ , shape :: {-# UNPACK #-} !Double+ } deriving (Eq)++instance Show WeibullDistribution where+ show (Weibull loc sc sh) = show "Weibull Distribution; loc: " ++ show loc ++ ", scale: " ++ show sc ++ " and shape: " ++ show sh++-- error message when initiating Weibull with scale and/or shape parameter less than 0.+weibullErrMsg :: Double -> Double -> Double -> String+weibullErrMsg loc scale sh = "Gev.Weibull: " + ++ "loc = " ++ show loc+ ++ " scale = " ++ show scale+ ++ " schape = " ++ show sh+ ++ ", but both the scale and shape parameters must be positive!"++weibullDistMaybe :: Double -> Double -> Double -> Maybe WeibullDistribution+weibullDistMaybe loc sc sh+ | sc > 0 && sh > 0 = Just $ Weibull loc sc sh+ | otherwise = Nothing++weibullDist :: Double -> Double -> Double -> WeibullDistribution+weibullDist loc sc sh = maybe (error $ weibullErrMsg loc sc sh) id $ weibullDistMaybe loc sh sc++-- | The CDF of the Weibull distribution+cdfWeibull :: WeibullDistribution -> Double -> Double+cdfWeibull (Weibull loc sc sh) x+ | x <= 0 = 0+ | otherwise = exp $ - ((- y) ** sh)+ where + y = (x - loc) / sc++-- | The PDF of the Weibull distribution+pdfWeibull :: WeibullDistribution -> Double -> Double+pdfWeibull (Weibull loc sc sh) x =+ let y = (x - loc) / sc+ const = sh / sc+ expterm = exp $ - ((- y) ** sh)+ middle = (- y) ** (sh - 1)+ in const * middle * expterm++quantileWeibull :: WeibullDistribution -> Double -> Double+quantileWeibull (Weibull loc sc sh) x+ | x == 0 = -inf+ | x == 1 = inf+ | x > 0 && x < 1 = loc - sc * (logexp ** pow)+ | otherwise =+ error $ "Gev.WeibullDistribution.quantile: The given value must be between 0 and 1, got: " ++ show x+ where+ inf = 1 / 0+ logexp = - log x+ pow = 1 / sh++instance Gev.Distribution WeibullDistribution where+ pdf = pdfWeibull+ cdf = cdfWeibull+ quantile = quantileWeibull
+ src/Tests/tests.hs view
@@ -0,0 +1,119 @@+-- | tests for the library++module Main where++import Gev+import Gev.Gumbel+import Gev.Frechet+import Gev.Weibull+import Gev.GevDist++import Test.HUnit+import qualified System.Exit as Exit++------------------------------------ Gumbel ------------------------------------++-- | Gumbel dist instance.+gumbel :: GumbelDistribution+gumbel = gumbelDist 0.5 2.0++-- | individual test+testGumbelCDF :: Test+testGumbelCDF = TestCase (assertEqual "Gumbel CDF test" 0.6235249162568004 (cdf gumbel 2.0))++testGumbelPDF :: Test+testGumbelPDF = TestCase (assertEqual "Gumbel PDF test" 0.14726615762017733 (pdf gumbel 2.0))++testGumbelQuant :: Test+testGumbelQuant = TestCase (assertEqual "Gumbel Qauntile test" 2.5618608663174456 (quantile gumbel 0.7))++------------------------------------ Frechet -----------------------------------++-- | Frechet dist instance.+frechet :: FrechetDistribution+frechet = frechetDist 1.0 0.1 1.0++-- | individual test+testFrechetCDF :: Test+testFrechetCDF = TestCase (assertEqual "Frechet CDF test" 0.951229424500714 (cdf frechet 3.0))++testFrechetPDF :: Test+testFrechetPDF = TestCase (assertEqual "Frechet PDF test" 0.023780735612517853 (pdf frechet 3.0))++testFrechetQuant :: Test+testFrechetQuant = TestCase (assertEqual "Frechet Qauntile test" 1.2803673252057128 (quantile frechet 0.7))++------------------------------------ Weibull -----------------------------------++-- | Weibull dist instance.+weibull :: WeibullDistribution+weibull = weibullDist 2.0 2.0 2.0++-- | individual test+testWeibullCDF :: Test+testWeibullCDF = TestCase (assertEqual "Weibull CDF test" 0.7788007830714049 (cdf weibull 1.0))++testWeibullPDF :: Test+testWeibullPDF = TestCase (assertEqual "Weibull PDF test" 0.38940039153570244 (pdf weibull 1.0))++testWeibullQuant :: Test+testWeibullQuant = TestCase (assertEqual "Weibull Qauntile test" 0.8055546158342233 (quantile weibull 0.7))+++------------------------------------ GEV -----------------------------------++--------- GEV with shape not equal to 0++-- | GEV dist instance.+gev :: GevDistribution+gev = gevDist 2.0 2.0 2.0++-- | individual test+testGEVCDF :: Test+testGEVCDF = TestCase (assertEqual "GEV CDF test" 0.4930686913952398 (cdf gev 3.0))++testGEVPDF :: Test+testGEVPDF = TestCase (assertEqual "GEV PDF test" 0.0871630538190878 (pdf gev 3.0))++testGEVQuant :: Test+testGEVQuant = TestCase (assertEqual "GEV Qauntile test" 8.860583704300595 (quantile gev 0.7))+++--------- GEV with shape equal to 0++-- | GEV dist instance with shape equal to 0.+gevAlt :: GevDistribution+gevAlt = gevDist 2.0 2.0 0.0++-- | individual test+testGEVCDFAlt :: Test+testGEVCDFAlt = TestCase (assertEqual "GEV CDF test" 0.545239211892605 (cdf gevAlt 3.0))++testGEVPDFAlt :: Test+testGEVPDFAlt = TestCase (assertEqual "GEV PDF test" 0.16535214944520904 (pdf gevAlt 3.0))++testGEVQuantAlt :: Test+testGEVQuantAlt = TestCase (assertEqual "GEV Qauntile test" 4.061860866317446 (quantile gevAlt 0.7))++-- | list of tests+testList :: Test+testList = TestList [ TestLabel "testGumbelCDF" testGumbelCDF+ , TestLabel "testGumbelPDF" testGumbelPDF+ , TestLabel "testGumbelQuant" testGumbelQuant+ , TestLabel "testFrechetCDF" testFrechetCDF+ , TestLabel "testFrechetPDF" testFrechetPDF+ , TestLabel "testFrechetQuant" testFrechetQuant+ , TestLabel "testWeibullCDF" testWeibullCDF+ , TestLabel "testWeibullPDF" testWeibullPDF+ , TestLabel "testWeibullQuant" testWeibullQuant+ , TestLabel "testGEVCDF" testGEVCDF+ , TestLabel "testGEVPDF" testGEVPDF+ , TestLabel "testGEVQuant" testGEVQuant+ , TestLabel "testGEVCDFAlt" testGEVCDFAlt+ , TestLabel "testGEVPDFAlt" testGEVPDFAlt+ , TestLabel "testGEVQuantAlt" testGEVQuantAlt]++main :: IO ()+main = do+ result <- runTestTT testList+ if failures result > 0 then Exit.exitFailure else Exit.exitSuccess