packages feed

hypergeometric 0.1.4.0 → 0.1.5.0

raw patch · 4 files changed

+56/−15 lines, 4 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

+ Math.Hypergeometric: euler :: (Ord a, Floating a) => a -> a -> a -> a -> a
+ Math.SpecialFunction: regbeta :: (Floating a, Ord a) => a -> a -> a -> a
- Math.Hypergeometric: erf :: (Eq a, Floating a) => a -> a
+ Math.Hypergeometric: erf :: (Ord a, Floating a) => a -> a
- Math.Hypergeometric: hypergeometric :: (Eq a, Fractional a) => [a] -> [a] -> a -> a
+ Math.Hypergeometric: hypergeometric :: (Ord a, Fractional a) => [a] -> [a] -> a -> a
- Math.Hypergeometric: ncdf :: (Eq a, Floating a) => a -> a
+ Math.Hypergeometric: ncdf :: (Ord a, Floating a) => a -> a

Files

CHANGELOG.md view
@@ -1,5 +1,11 @@ # hypergeometric +## 0.1.5.0++  * Throw errors when outside of the radius of convergence for a series+  * Add `euler` which transforms ₂F₁ to something that converges for all z.+  * Fix egregious typo in documentation+ ## 0.1.4.0    * Add `agm`, `completeElliptic`
hypergeometric.cabal view
@@ -1,6 +1,6 @@ cabal-version:   1.18 name:            hypergeometric-version:         0.1.4.0+version:         0.1.5.0 license:         AGPL-3 license-file:    COPYING copyright:       Copyright: (c) 2022 Vanessa McHale@@ -9,7 +9,7 @@ bug-reports:     https://github.com/vmchale/hypergeometric/issues synopsis:        Hypergeometric functions description:-    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.+    Haskell implementation of hypergeometric functions and associated statistical and special functions: erf, normal cdf, incomplete beta, regularized beta, F-distribution cdf, \(\chi^2\)-distribution cdf, t-distrubtion cdf.     Also includes Lanczos' approximation of the gamma function.  category:        Math, Statistics
src/Math/Hypergeometric.hs view
@@ -1,5 +1,6 @@--- | See McHale, Vanessa [\"Hypergeometric Functions for Statistical Computing\"](http://vmchale.com/static/serve/hypergeometric.pdf) and especially Shaw, Ernest [\"Hypergeometric Functions and CDFs in J\"](https://www.jsoftware.com/papers/jhyper.pdf)+-- | See Shaw, Ewart [\"Hypergeometric Functions and CDFs in J\"](https://www.jsoftware.com/papers/jhyper.pdf). module Math.Hypergeometric ( hypergeometric+                           , euler                            , erf                            , ncdf                            ) where@@ -16,27 +17,56 @@ {-# SPECIALIZE ncdf :: Double -> Double #-} {-# SPECIALIZE ncdf :: Float -> Float #-} -- | CDF of the standard normal \( N(0,1) \)-ncdf :: (Eq a, Floating a) => a -> a+ncdf :: (Ord a, Floating a) => a -> a ncdf z = (1/2) * (1 + erf (z / sqrt 2))  {-# SPECIALIZE erf :: Double -> Double #-} {-# SPECIALIZE erf :: Float -> Float #-} -- | [erf](https://mathworld.wolfram.com/Erf.html)-erf :: (Eq a, Floating a) => a -> a-erf z = (2 * z * exp (-z^(2::Int)) / sqrt pi) * hypergeometric [1] [3/2] (z^(2::Int))+erf :: (Ord a, Floating a) => a -> a+erf z = (2 * z * exp (-(z^(2::Int))) / sqrt pi) * hypergeometric [1] [3/2] (z^(2::Int)) +{-# SPECIALIZE euler :: Double -> Double -> Double -> Double -> Double #-}+{-# SPECIALIZE euler :: Float -> Float -> Float -> Float -> Float #-}+-- | Euler's transform:+--+-- \( \displaystyle _2F_1(a,b;c;z) = (1-z)^{-a} {}_2F_1\left(a,c-b;c;\frac{z}{z-1}\right) \)+--+-- The right-hand side converges for all \(z\).+--+-- Koekoek, Roelef and Swarttouw, René F. [The Askey-scheme of hypergeometric orthogonal polynomials and its q-analogue](https://arxiv.org/abs/math/9602214).+--+-- @since 0.1.5.0+euler :: (Ord a, Floating a)+      => a -- ^ \(a\)+      -> a -- ^ \(b\)+      -> a -- ^ \(c\)+      -> a -- ^ \(z\)+      -> a+euler a b c z = hypergeometric [a, c-b] [c] (z/z-1)+ {-# SPECIALIZE hypergeometric :: [Double] -> [Double] -> Double -> Double #-} {-# SPECIALIZE hypergeometric :: [Float] -> [Float] -> Float -> Float #-} -- | \( _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!} \) --+-- The radius of convergence is+--+-- \( \rho = \begin{cases} \infty & \text{if} & p<q+1 \\ 1 & \text{if} & p=q+1 \\ 0 & \text{if} & p>q+1 \\ \end{cases} \)+-- -- This iterates until the result stabilizes.-hypergeometric :: (Eq a, Fractional a)+hypergeometric :: (Ord a, Fractional a)                => [a] -- ^ \( a_1,\ldots,a_p \)                -> [a] -- ^ \( b_1,\ldots,b_q \)                -> a -- ^ \( z \)                -> a hypergeometric as bs z = sumUntilEq-    [ (product (fmap (`risingFactorial` n) as) / product (fmap (`risingFactorial` n) bs)) * (z ^ n) / factorial n | n <- [0..] ]+    [ (product (fmap (`risingFactorial` n) as) / product (fmap (`risingFactorial` n) bs)) * (zϵ ^ n) / factorial n | n <- [0..] ]+  where+    zϵ = 𝜌 z+    𝜌 = case compare (length as) (length bs+1) of+        LT -> id+        EQ -> if z > 1 then error "Outside the radius of convergence." else id+        GT -> if z /= 0 then error "Outside the radius of convergence." else id  sumUntilEq :: (Eq a, Num a) => [a] -> a sumUntilEq = sumUntilEqLoop 0
src/Math/SpecialFunction.hs view
@@ -1,5 +1,6 @@ module Math.SpecialFunction ( incbeta                             , beta+                            , regbeta                             , bessel1                             , gamma                             , gammaln@@ -14,14 +15,14 @@  {-# SPECIALIZE tcdf :: Double -> Double -> Double #-} {-# SPECIALIZE tcdf :: Float -> Float -> Float #-}--- | Converges if and only if \(|x| < \sqrt{\nu} \)+-- | Converges if and only if \(|x| \leq \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)/𝜈)+tcdf 𝜈 x = 0.5 + x * gamma (0.5*(𝜈+1)) / (sqrt(pi*𝜈) * gamma(𝜈/2)) * hypergeometric [0.5, 0.5*(𝜈+1)] [1.5] (-(x^(2::Int))/𝜈)  {-# SPECIALIZE bessel1 :: Double -> Double -> Double #-} {-# SPECIALIZE bessel1 :: Float -> Float -> Float #-}@@ -69,15 +70,13 @@ -- | \(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 #-} {-# SPECIALIZE incbeta :: Float -> Float -> Float -> Float #-}--- | Incomplete beta function, \(|z|<1\)+-- | Incomplete beta function, \(|z|\leq 1\) -- -- Calculated with \(B(z;a,b)=\displaystyle\frac{z^a}{a}{}_2F_1(a, 1-b; a+1; z)\) --@@ -91,8 +90,14 @@  {-# SPECIALIZE regbeta :: Double -> Double -> Double -> Double #-} {-# SPECIALIZE regbeta :: Float -> Float -> Float -> Float #-}--- | \(I(z;a,b) = \displaystyle\frac{B(z;a,b)}{B(a,b)}\)-regbeta :: (Floating a, Ord a) => a -> a -> a -> a+-- | [Regularized beta function](https://mathworld.wolfram.com/RegularizedBetaFunction.html), \(|z|\leq 1\)+--+-- \(I(z;a,b) = \displaystyle\frac{B(z;a,b)}{B(a,b)}\)+regbeta :: (Floating a, Ord a) +        => a -- ^ \(z\)+        -> a -- ^ \(a\)+        -> a -- ^ \(b\)+        -> a regbeta z a b = incbeta z a b / beta a b  {-# SPECIALIZE fcdf :: Double -> Double -> Double -> Double #-}