hypergeometric 0.1.1.0 → 0.1.2.0

raw patch · 4 files changed

+57/−19 lines, 4 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• CHANGELOG.md +4/−0
• hypergeometric.cabal +3/−2
• src/Math/Hypergeometric.hs +4/−8
• src/Math/SpecialFunction.hs +46/−9

@@ -1,5 +1,9 @@ # hypergeometric +## 0.1.2.0++  * Add `tcdf`, `chisqcdf`, `fcdf`+ ## 0.1.1.0    * Add `Math.SpecialFunction` with `gamma`, `beta` etc.

@@ -1,6 +1,6 @@ cabal-version:   1.18 name:            hypergeometric-version:         0.1.1.0+version:         0.1.2.0 license:         AGPL-3 license-file:    COPYING copyright:       Copyright: (c) 2022 Vanessa McHale@@ -9,7 +9,8 @@ bug-reports:     https://github.com/vmchale/hypergeometric/issues synopsis:        Hypergeometric functions description:-    Haskell implementation of hypergeometric functions and associated statistical functions, viz. erf, normal cdf+    Haskell implementation of hypergeometric functions and associated statistical and special functions, viz. erf, normal cdf, incomplete beta, F-distribution cdf, \(\chi^2\)-distribution cdf, t-distrubtion cdf.+    Also includes Lanczos' approximation of the gamma function.  category:        Math, Statistics build-type:      Simple

@@ -6,9 +6,6 @@  import           Data.Functor ((<$>)) --- choose :: Integral a => a -> a -> a--- choose n k = product [(n-k+1) .. n] `quot` factorial (fromIntegral k)- risingFactorial :: Num a => a -> Int -> a risingFactorial _ 0 = 1 risingFactorial a n = (a + fromIntegral n - 1) * risingFactorial a (n-1)@@ -16,6 +13,9 @@ factorial :: Num a => Int -> a factorial n = product (fromIntegral <$> [1..n]) +-- prop_cdf :: (Double -> Double) -> Double -> Bool+-- prop_cdf f x = f x <= 1+ {-# SPECIALIZE ncdf :: Double -> Double #-} -- | CDF of the standard normal \( N(0,1) \) ncdf :: (Eq a, Floating a) => a -> a@@ -29,8 +29,7 @@ {-# SPECIALIZE hypergeometric :: [Double] -> [Double] -> Double -> Double #-} -- | \( _pF_q(a_1,\ldots,a_p;b_1,\ldots,b_q;z) = \displaystyle\sum_{n=0}^\infty\frac{(a_1)_n\cdots(a_p)_n}{(b_1)_b\cdots(b_q)_n}\frac{z^n}{n!} \) ----- This iterates until the result stabilizes, so don't use it on--- arbitrary-precision types!+-- This iterates until the result stabilizes. hypergeometric :: (Eq a, Fractional a)                => [a] -- ^ \( a_1,\ldots,a_p \)                -> [a] -- ^ \( b_1,\ldots,b_q \)@@ -38,9 +37,6 @@                -> a hypergeometric as bs z = sumUntilEq     [ (product (fmap (`risingFactorial` n) as) / product (fmap (`risingFactorial` n) bs)) * (z ^ n) / factorial n | n <- [0..] ]-    -- [ exp (nth n) | n <- [0..] ]-    -- where nth n = sum (fmap (log . (`risingFactorial` n)) as) - sum (fmap (log . (`risingFactorial` n)) bs) + fromIntegral n * log z - log (factorial n)-    -- TODO: Revisit the exponential approach using complex numbers?  sumUntilEq :: (Eq a, Num a) => [a] -> a sumUntilEq = sumUntilEqLoop 0

@@ -2,26 +2,68 @@                             , beta                             , gamma                             , gammaln+                            , fcdf+                            , chisqcdf+                            , tcdf                             ) where  import           Math.Hypergeometric --- prop_betamatch :: Double -> Double -> Bool--- prop_betamatch x y = x <= 0 || y <= 0 || abs (beta x y - incbeta 1 x y) < 1e-15+{-# SPECIALIZE tcdf :: Double -> Double -> Double #-}+-- | Converges if and only if \(|x| < \sqrt{\nu} \)+--+-- @since 0.1.2.0+tcdf :: (Floating a, Ord a)+     => a -- ^ \(\nu\) (degrees of freedom)+     -> a -- ^ \(x\)+     -> a+tcdf 𝜈 x = 0.5 + x * gamma (0.5*(𝜈+1)) / (sqrt(pi*𝜈) * gamma(𝜈/2)) * hypergeometric [0.5, 0.5*(𝜈+1)] [1.5] (-x^(2::Int)/𝜈) +-- | @since 0.1.2.0+{-# SPECIALIZE chisqcdf :: Double -> Double -> Double #-}+chisqcdf :: (Floating a, Ord a)+         => a -- ^ \(r\) (degrees of freedom)+         -> a -- ^ \(\chi^2\)+         -> a+chisqcdf r x = incgamma (0.5*r) (0.5*x) / gamma (0.5*r)++{-# SPECIALIZE incgamma :: Double -> Double -> Double #-}+-- | \(a^{-1}x^a{}_1F_1(a;1+a;-x) \)+incgamma :: (Floating a, Ord a) => a -> a -> a+incgamma a x = (1/a) * x ** a * hypergeometric [a] [1+a] (-x)+-- TODO: writeup?+--+-- chisqcdf 10 28 works better this way than w/ e^-x ... x+--+-- (1 2 H. _1.1) 1 hangs indefinitely+ {-# SPECIALIZE incbeta :: Double -> Double -> Double -> Double #-}--- | Incomplete beta function.+-- | Incomplete beta function, \(|z|<1\) -- -- Calculated with \(B(z;a,b)=\displaystyle\frac{z^a}{a}{}_2F_1(a, 1-b; a+1; z)\) -- -- @since 0.1.1.0-incbeta :: (Floating a, Eq a)+incbeta :: (Floating a, Ord a)         => a -- ^ \(z\)         -> a -- ^ \(a\)         -> a -- ^ \(b\)         -> a incbeta z a b = z**a/a * hypergeometric [a,1-b] [a+1] z +{-# SPECIALIZE regbeta :: Double -> Double -> Double -> Double #-}+-- | \(I(z;a,b) = \displaystyle\frac{B(z;a,b)}{B(a,b)}\)+regbeta :: (Floating a, Ord a) => a -> a -> a -> a+regbeta z a b = incbeta z a b / beta a b++{-# SPECIALIZE fcdf :: Double -> Double -> Double -> Double #-}+-- | @since 0.1.2.0+fcdf :: (Floating a, Ord a)+     => a -- ^ \(n\)+     -> a -- ^ \(m\)+     -> a -- ^ \(x\)+     -> a+fcdf n m x = regbeta (n * x / (m + n * x)) (0.5 * n) (0.5 * m) -- we can use hypergeo because nx/(m+nx) < 1+ {-# SPECIALIZE beta :: Double -> Double -> Double #-} -- | \(B(x, y) = \displaystyle\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}\) --@@ -40,11 +82,6 @@ -- @since 0.1.1.0 gamma :: (Floating a, Ord a) => a -> a gamma = exp . gammaln---- gamma from beta:--- Γ(z)Γ(1-z) = 𝜋/sin(𝜋z)------ THENCE, B(z,1-z)=Γ(z)Γ(1-z)/Γ(1)=...  {-# SPECIALIZE gammaln :: Double -> Double #-} -- | \(\text{log} (\Gamma(z))\)