hypergeometric 0.1.0.0 → 0.1.1.0

raw patch · 4 files changed

+100/−2 lines, 4 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• CHANGELOG.md +4/−0
• hypergeometric.cabal +6/−2
• src/Math/Hypergeometric.hs +9/−0
• src/Math/SpecialFunction.hs +81/−0

@@ -1,5 +1,9 @@ # hypergeometric +## 0.1.1.0++  * Add `Math.SpecialFunction` with `gamma`, `beta` etc.+ ## 0.1.0.0  Initial release

@@ -1,11 +1,12 @@ cabal-version:   1.18 name:            hypergeometric-version:         0.1.0.0+version:         0.1.1.0 license:         AGPL-3 license-file:    COPYING copyright:       Copyright: (c) 2022 Vanessa McHale maintainer:      vamchale@gmail.com author:          Vanessa McHale+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@@ -21,7 +22,10 @@     location: https://github.com/vmchale/hypergeometric  library-    exposed-modules:  Math.Hypergeometric+    exposed-modules:+        Math.Hypergeometric+        Math.SpecialFunction+     hs-source-dirs:   src     default-language: Haskell2010     ghc-options:      -Wall

@@ -6,6 +6,9 @@  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)@@ -25,6 +28,9 @@  {-# 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! hypergeometric :: (Eq a, Fractional a)                => [a] -- ^ \( a_1,\ldots,a_p \)                -> [a] -- ^ \( b_1,\ldots,b_q \)@@ -32,6 +38,9 @@                -> 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

@@ -0,0 +1,81 @@+module Math.SpecialFunction ( incbeta+                            , beta+                            , gamma+                            , gammaln+                            ) 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 incbeta :: Double -> Double -> Double -> Double #-}+-- | Incomplete beta function.+--+-- 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)+        => a -- ^ \(z\)+        -> a -- ^ \(a\)+        -> a -- ^ \(b\)+        -> a+incbeta z a b = z**a/a * hypergeometric [a,1-b] [a+1] z++{-# SPECIALIZE beta :: Double -> Double -> Double #-}+-- | \(B(x, y) = \displaystyle\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}\)+--+-- This uses 'gammaln' under the hood to extend its domain somewhat.+--+-- @since 0.1.1.0+beta :: (Floating a, Ord a) => a -> a -> a+beta x y = exp (betaln x y)++{-# SPECIALIZE betaln :: Double -> Double -> Double #-}+betaln :: (Floating a, Ord a) => a -> a -> a+betaln x y = gammaln x + gammaln y - gammaln (x+y)++-- | \(\Gamma(z)\)+--+-- @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))\)+--+-- Lanczos approximation.+-- This is exactly the approach described in Press, William H. et al. /Numerical Recipes/, 3rd ed., extended to work on negative real numbers.+--+-- @since 0.1.1.0+gammaln :: (Floating a, Ord a)+        => a -- ^ \( z \)+        -> a+gammaln 0 = -log 0+gammaln z | z >= 0.5 = (z' + 1/2) * log (z' + 𝛾 + 1/2) - (z' + 1/2 + 𝛾) + log (sqrt (2*pi) * (c0 + sum series))+    where series = zipWith (\c x -> c / (z' + fromIntegral x)) coeff [(1::Int)..]+          c0 = 0.999999999999997092+          -- constants from Numerical Recipes+          coeff = [ 57.1562356658629235+                  , -59.5979603554754912+                  , 14.1360979747417471+                  , -0.491913816097620199+                  , 0.339946499848118887e-4+                  , 0.465236289270485756e-4+                  , -0.983744753048795646e-4+                  , 0.158088703224912494e-3+                  , -0.210264441724104883e-3+                  , 0.217439618115212643e-3+                  , -0.164318106536763890e-3+                  , 0.844182239838527433e-4+                  , -0.261908384015814087e-4+                  , 0.368991826595316234e-5+                  ]+          𝛾 = 607/128+          z' = z-1+gammaln z = log pi - log (sin (pi * z)) - gammaln (1 - z)