math-functions 0.1.7.0 → 0.2.0.0
raw patch · 13 files changed
+969/−477 lines, 13 filesdep ~basedep ~erf
Dependency ranges changed: base, erf
Files
- Numeric/MathFunctions/Comparison.hs +32/−5
- Numeric/MathFunctions/Constants.hs +13/−1
- Numeric/Polynomial/Chebyshev.hs +6/−3
- Numeric/RootFinding.hs +173/−0
- Numeric/Series.hs +180/−0
- Numeric/SpecFunctions.hs +11/−4
- Numeric/SpecFunctions/Extra.hs +62/−2
- Numeric/SpecFunctions/Internal.hs +420/−238
- Numeric/Sum.hs +1/−1
- changelog.md +33/−0
- math-functions.cabal +20/−5
- tests/Tests/SpecFunctions.hs +18/−12
- tests/Tests/SpecFunctions_flymake.hs +0/−206
Numeric/MathFunctions/Comparison.hs view
@@ -19,6 +19,7 @@ -- * Ulps-based comparison , addUlps , ulpDistance+ , ulpDelta , within ) where @@ -36,7 +37,7 @@ -- | -- Calculate relative error of two numbers: ----- > |a - b| / max |a| |b|+-- \[ \frac{|a - b|}{\max(|a|,|b|)} \] -- -- It lies in [0,1) interval for numbers with same sign and (1,2] for -- numbers with different sign. If both arguments are zero or negative@@ -49,9 +50,9 @@ -- | Check that relative error between two numbers @a@ and @b@. If -- 'relativeError' returns NaN it returns @False@.-eqRelErr :: Double -- ^ @eps@ relative error should be in [0,1) range- -> Double -- ^ @a@- -> Double -- ^ @b@+eqRelErr :: Double -- ^ /eps/ relative error should be in [0,1) range+ -> Double -- ^ /a/+ -> Double -- ^ /b/ -> Bool eqRelErr eps a b = relativeError a b < eps @@ -80,7 +81,8 @@ -- | -- Measure distance between two @Double@s in ULPs (units of least--- precision).+-- precision). Note that it's different from @abs (ulpDelta a b)@+-- since it returns correct result even when 'ulpDelta' overflows. ulpDistance :: Double -> Double -> Word64@@ -100,6 +102,31 @@ d | ai > bi = ai - bi | otherwise = bi - ai return $! d++-- |+-- Measure signed distance between two @Double@s in ULPs (units of least+-- precision). Note that unlike 'ulpDistance' it can overflow.+--+-- > >>> ulpDelta 1 (1 + m_epsilon)+-- > 1+ulpDelta :: Double+ -> Double+ -> Int64+ulpDelta a b = runST $ do+ buf <- newByteArray 8+ ai0 <- writeByteArray buf 0 a >> readByteArray buf 0+ bi0 <- writeByteArray buf 0 b >> readByteArray buf 0+ -- IEEE754 floats use most significant bit as sign bit (not+ -- 2-complement) and we need to rearrange representations of float+ -- number so that they could be compared lexicographically as+ -- Word64.+ let big = 0x8000000000000000 :: Word64+ order i | i < big = i + big+ | otherwise = maxBound - i+ ai = order ai0+ bi = order bi0+ return $! fromIntegral $ bi - ai+ -- | Compare two 'Double' values for approximate equality, using -- Dawson's method.
Numeric/MathFunctions/Constants.hs view
@@ -19,6 +19,8 @@ , m_pos_inf , m_neg_inf , m_NaN+ , m_max_log+ , m_min_log -- * Mathematical constants , m_1_sqrt_2 , m_2_sqrt_pi@@ -32,11 +34,12 @@ -- IEE754 constants ---------------------------------------------------------------- --- | A very large number.+-- | Largest representable finite value. m_huge :: Double m_huge = 1.7976931348623157e308 {-# INLINE m_huge #-} +-- | The smallest representable positive normalized value. m_tiny :: Double m_tiny = 2.2250738585072014e-308 {-# INLINE m_tiny #-}@@ -61,6 +64,15 @@ m_NaN = 0/0 {-# INLINE m_NaN #-} +-- | Maximum possible finite value of @log x@+m_max_log :: Double+m_max_log = 709.782712893384+{-# INLINE m_max_log #-}++-- | Logarithm of smallest normalized double ('m_tiny')+m_min_log :: Double+m_min_log = -708.3964185322641+{-# INLINE m_min_log #-} ----------------------------------------------------------------
Numeric/Polynomial/Chebyshev.hs view
@@ -27,9 +27,12 @@ -- A Chebyshev polynomial of the first kind is defined by the -- following recurrence: ----- > t 0 _ = 1--- > t 1 x = x--- > t n x = 2 * x * t (n-1) x - t (n-2) x+-- \[\begin{aligned}+-- T_0(x) &= 1 \\+-- T_1(x) &= x \\+-- T_{n+1}(x) &= 2xT_n(x) - T_{n-1}(x) \\+-- \end{aligned}+-- \] data C = C {-# UNPACK #-} !Double {-# UNPACK #-} !Double
+ Numeric/RootFinding.hs view
@@ -0,0 +1,173 @@+{-# LANGUAGE BangPatterns, DeriveDataTypeable, DeriveGeneric, CPP #-}+-- |+-- Module : Numeric.RootFinding+-- Copyright : (c) 2011 Bryan O'Sullivan+-- License : BSD3+--+-- Maintainer : bos@serpentine.com+-- Stability : experimental+-- Portability : portable+--+-- Haskell functions for finding the roots of real functions of real arguments.+module Numeric.RootFinding+ (+ Root(..)+ , fromRoot+ , ridders+ , newtonRaphson+ -- * References+ -- $references+ ) where++import Control.Applicative (Alternative(..), Applicative(..))+import Control.Monad (MonadPlus(..), ap)+import Data.Data (Data, Typeable)+#if __GLASGOW_HASKELL__ > 704+import GHC.Generics (Generic)+#endif+import Numeric.MathFunctions.Comparison (within)+++-- | The result of searching for a root of a mathematical function.+data Root a = NotBracketed+ -- ^ The function does not have opposite signs when+ -- evaluated at the lower and upper bounds of the search.+ | SearchFailed+ -- ^ The search failed to converge to within the given+ -- error tolerance after the given number of iterations.+ | Root a+ -- ^ A root was successfully found.+ deriving (Eq, Read, Show, Typeable, Data+#if __GLASGOW_HASKELL__ > 704+ , Generic+#endif+ )+++instance Functor Root where+ fmap _ NotBracketed = NotBracketed+ fmap _ SearchFailed = SearchFailed+ fmap f (Root a) = Root (f a)++instance Monad Root where+ NotBracketed >>= _ = NotBracketed+ SearchFailed >>= _ = SearchFailed+ Root a >>= m = m a++ return = Root++instance MonadPlus Root where+ mzero = SearchFailed++ r@(Root _) `mplus` _ = r+ _ `mplus` p = p++instance Applicative Root where+ pure = Root+ (<*>) = ap++instance Alternative Root where+ empty = SearchFailed++ r@(Root _) <|> _ = r+ _ <|> p = p++-- | Returns either the result of a search for a root, or the default+-- value if the search failed.+fromRoot :: a -- ^ Default value.+ -> Root a -- ^ Result of search for a root.+ -> a+fromRoot _ (Root a) = a+fromRoot a _ = a+++-- | Use the method of Ridders to compute a root of a function.+--+-- The function must have opposite signs when evaluated at the lower+-- and upper bounds of the search (i.e. the root must be bracketed).+ridders :: Double -- ^ Absolute error tolerance.+ -> (Double,Double) -- ^ Lower and upper bounds for the search.+ -> (Double -> Double) -- ^ Function to find the roots of.+ -> Root Double+ridders tol (lo,hi) f+ | flo == 0 = Root lo+ | fhi == 0 = Root hi+ | flo*fhi > 0 = NotBracketed -- root is not bracketed+ | otherwise = go lo flo hi fhi 0+ where+ go !a !fa !b !fb !i+ -- Root is bracketed within 1 ulp. No improvement could be made+ | within 1 a b = Root a+ -- Root is found. Check that f(m) == 0 is nessesary to ensure+ -- that root is never passed to 'go'+ | fm == 0 = Root m+ | fn == 0 = Root n+ | d < tol = Root n+ -- Too many iterations performed. Fail+ | i >= (100 :: Int) = SearchFailed+ -- Ridder's approximation coincide with one of old+ -- bounds. Revert to bisection+ | n == a || n == b = case () of+ _| fm*fa < 0 -> go a fa m fm (i+1)+ | otherwise -> go m fm b fb (i+1)+ -- Proceed as usual+ | fn*fm < 0 = go n fn m fm (i+1)+ | fn*fa < 0 = go a fa n fn (i+1)+ | otherwise = go n fn b fb (i+1)+ where+ d = abs (b - a)+ dm = (b - a) * 0.5+ !m = a + dm+ !fm = f m+ !dn = signum (fb - fa) * dm * fm / sqrt(fm*fm - fa*fb)+ !n = m - signum dn * min (abs dn) (abs dm - 0.5 * tol)+ !fn = f n+ !flo = f lo+ !fhi = f hi+++-- | Solve equation using Newton-Raphson iterations.+--+-- This method require both initial guess and bounds for root. If+-- Newton step takes us out of bounds on root function reverts to+-- bisection.+newtonRaphson+ :: Double+ -- ^ Required precision+ -> (Double,Double,Double)+ -- ^ (lower bound, initial guess, upper bound). Iterations will no+ -- go outside of the interval+ -> (Double -> (Double,Double))+ -- ^ Function to finds roots. It returns pair of function value and+ -- its derivative+ -> Root Double+newtonRaphson !prec (!low,!guess,!hi) function+ = go low guess hi+ where+ go !xMin !x !xMax+ | f == 0 = Root x+ | abs (dx / x) < prec = Root x+ | otherwise = go xMin' x' xMax'+ where+ (f,f') = function x+ -- Calculate Newton-Raphson step+ delta | f' == 0 = error "handle f'==0"+ | otherwise = f / f'+ -- Calculate new approximation and actual change of approximation+ (dx,x') | z <= xMin = let d = 0.5*(x - xMin) in (d, x - d)+ | z >= xMax = let d = 0.5*(x - xMax) in (d, x - d)+ | otherwise = (delta, z)+ where z = x - delta+ -- Update root bracket+ xMin' | dx < 0 = x+ | otherwise = xMin+ xMax' | dx > 0 = x+ | otherwise = xMax++++-- $references+--+-- * Ridders, C.F.J. (1979) A new algorithm for computing a single+-- root of a real continuous function.+-- /IEEE Transactions on Circuits and Systems/ 26:979–980.
+ Numeric/Series.hs view
@@ -0,0 +1,180 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE ExistentialQuantification #-}+-- |+-- Module : Numeric.Series+-- Copyright : (c) 2016 Alexey Khudyakov+-- License : BSD3+--+-- Maintainer : alexey.skladnoy@gmail.com, bos@serpentine.com+-- Stability : experimental+-- Portability : portable+--+-- Functions for working with infinite sequences. In particular+-- summation of series and evaluation of continued fractions.+module Numeric.Series (+ -- * Data type for infinite sequences.+ Sequence(..)+ -- * Constructors+ , enumSequenceFrom+ , enumSequenceFromStep+ , scanSequence+ -- * Summation of series+ , sumSeries+ , sumPowerSeries+ , sequenceToList+ -- * Evaluation of continued fractions+ , evalContFractionB+ ) where++import Control.Applicative+import Data.List (unfoldr)++import Numeric.MathFunctions.Constants (m_epsilon)+++----------------------------------------------------------------++-- | Infinite series. It's represented as opaque state and step+-- function.+data Sequence a = forall s. Sequence s (s -> (a,s))++instance Functor Sequence where+ fmap f (Sequence s0 step) = Sequence s0 (\s -> let (a,s') = step s in (f a, s'))+ {-# INLINE fmap #-}++instance Applicative Sequence where+ pure a = Sequence () (\() -> (a,()))+ Sequence sA fA <*> Sequence sB fB = Sequence (sA,sB) $ \(!sa,!sb) ->+ let (a,sa') = fA sa+ (b,sb') = fB sb+ in (a b, (sa',sb'))+ {-# INLINE pure #-}+ {-# INLINE (<*>) #-}++-- | Elementwise operations with sequences+instance Num a => Num (Sequence a) where+ (+) = liftA2 (+)+ (*) = liftA2 (*)+ (-) = liftA2 (-)+ {-# INLINE (+) #-}+ {-# INLINE (*) #-}+ {-# INLINE (-) #-}+ abs = fmap abs+ signum = fmap signum+ fromInteger = pure . fromInteger+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ {-# INLINE fromInteger #-}++-- | Elementwise operations with sequences+instance Fractional a => Fractional (Sequence a) where+ (/) = liftA2 (/)+ recip = fmap recip+ fromRational = pure . fromRational+ {-# INLINE (/) #-}+ {-# INLINE recip #-}+ {-# INLINE fromRational #-}++++----------------------------------------------------------------+-- Constructors+----------------------------------------------------------------++-- | @enumSequenceFrom x@ generate sequence:+--+-- \[ a_n = x + n \]+enumSequenceFrom :: Num a => a -> Sequence a+enumSequenceFrom i = Sequence i (\n -> (n,n+1))+{-# INLINE enumSequenceFrom #-}++-- | @enumSequenceFromStep x d@ generate sequence:+--+-- \[ a_n = x + nd \]+enumSequenceFromStep :: Num a => a -> a -> Sequence a+enumSequenceFromStep n d = Sequence n (\i -> (i,i+d))+{-# INLINE enumSequenceFromStep #-}++-- | Analog of 'scanl' for sequence.+scanSequence :: (b -> a -> b) -> b -> Sequence a -> Sequence b+{-# INLINE scanSequence #-}+scanSequence f b0 (Sequence s0 step) = Sequence (b0,s0) $ \(b,s) ->+ let (a,s') = step s+ b' = f b a+ in (b,(b',s'))+++----------------------------------------------------------------+-- Evaluation of series+----------------------------------------------------------------++-- | Calculate sum of series+--+-- \[ \sum_{i=0}^\infty a_i \]+--+-- Summation is stopped when+--+-- \[ a_{n+1} < \varepsilon\sum_{i=0}^n a_i \]+--+-- where ε is machine precision ('m_epsilon')+sumSeries :: Sequence Double -> Double+{-# INLINE sumSeries #-}+sumSeries (Sequence sInit step)+ = go x0 s0+ where + (x0,s0) = step sInit+ go x s | abs (d/x) >= m_epsilon = go x' s'+ | otherwise = x'+ where (d,s') = step s+ x' = x + d++-- | Calculate sum of series+--+-- \[ \sum_{i=0}^\infty x^ia_i \]+--+-- Calculation is stopped when next value in series is less than+-- ε·sum.+sumPowerSeries :: Double -> Sequence Double -> Double+sumPowerSeries x ser = sumSeries $ liftA2 (*) (scanSequence (*) 1 (pure x)) ser+{-# INLINE sumPowerSeries #-}++-- | Convert series to infinite list+sequenceToList :: Sequence a -> [a]+sequenceToList (Sequence s f) = unfoldr (Just . f) s++++----------------------------------------------------------------+-- Evaluation of continued fractions+----------------------------------------------------------------++-- |+-- Evaluate continued fraction using modified Lentz algorithm.+-- Sequence contain pairs (a[i],b[i]) which form following expression:+--+-- \[+-- b_0 + \frac{a_1}{b_1+\frac{a_2}{b_2+\frac{a_3}{b_3 + \cdots}}}+-- \]+--+-- Modified Lentz algorithm is described in Numerical recipes 5.2+-- "Evaluation of Continued Fractions"+evalContFractionB :: Sequence (Double,Double) -> Double+{-# INLINE evalContFractionB #-}+evalContFractionB (Sequence sInit step)+ = let ((_,b0),s0) = step sInit+ f0 = maskZero b0+ in go f0 f0 0 s0+ where+ tiny = 1e-60+ maskZero 0 = tiny+ maskZero x = x+ + go f c d s+ | abs (delta - 1) >= m_epsilon = go f' c' d' s'+ | otherwise = f'+ where+ ((a,b),s') = step s+ d' = recip $ maskZero $ b + a*d+ c' = maskZero $ b + a/c + delta = c'*d'+ f' = f*delta
Numeric/SpecFunctions.hs view
@@ -1,4 +1,3 @@-{-# LANGUAGE BangPatterns, ScopedTypeVariables #-} -- | -- Module : Numeric.SpecFunctions -- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan@@ -30,7 +29,10 @@ , sinc -- * Logarithm , log1p+ , log1pmx , log2+ -- * Exponent+ , expm1 -- * Factorial , factorial , logFactorial@@ -81,6 +83,11 @@ -- Vol. 22, No. 3 (1973), pp. 411-414 -- <http://www.jstor.org/pss/2346798> ----- * Shea, B. (1988) Algorithm AS 239: Chi-squared and incomplete--- gamma integral. /Applied Statistics/--- 37(3):466–473. <http://www.jstor.org/stable/2347328>+-- * Temme, N.M. (1992) Asymptotic inversion of the incomplete beta+-- function. /Journal of Computational and Applied Mathematics+-- 41(1992) 145-157.+--+-- * Temme, N.M. (1994) A set of algorithms for the incomplete gamma+-- functions. /Probability in the Engineering and Informational+-- Sciences/, 8, 1994, 291-307. Printed in the U.S.A.+
Numeric/SpecFunctions/Extra.hs view
@@ -12,10 +12,12 @@ bd0 , chooseExact , logChooseFast+ , logGammaAS245+ , logGammaCorrection ) where -import Numeric.MathFunctions.Constants (m_NaN)-import Numeric.SpecFunctions.Internal (chooseExact,logChooseFast)+import Numeric.MathFunctions.Constants (m_NaN,m_pos_inf)+import Numeric.SpecFunctions.Internal (chooseExact,logChooseFast,logGammaCorrection) -- | Evaluate the deviance term @x log(x/np) + np - x@. bd0 :: Double -- ^ @x@@@ -34,3 +36,61 @@ loop j ej s = case s + ej/(2*j+1) of s' | s' == s -> s' -- FIXME: Comparing Doubles for equality! | otherwise -> loop (j+1) (ej*vv) s'++++-- | Compute the logarithm of the gamma function Γ(/x/). Uses+-- Algorithm AS 245 by Macleod.+--+-- Gives an accuracy of 10-12 significant decimal digits, except+-- for small regions around /x/ = 1 and /x/ = 2, where the function+-- goes to zero. For greater accuracy, use 'logGammaL'.+--+-- Returns ∞ if the input is outside of the range (0 < /x/ ≤ 1e305).+logGammaAS245 :: Double -> Double+-- Adapted from http://people.sc.fsu.edu/~burkardt/f_src/asa245/asa245.html+logGammaAS245 x+ | x <= 0 = m_pos_inf+ -- Handle positive infinity. logGamma overflows before 1e308 so+ -- it's safe+ | x > 1e308 = m_pos_inf+ -- Normal cases+ | x < 1.5 = a + c *+ ((((r1_4 * b + r1_3) * b + r1_2) * b + r1_1) * b + r1_0) /+ ((((b + r1_8) * b + r1_7) * b + r1_6) * b + r1_5)+ | x < 4 = (x - 2) *+ ((((r2_4 * x + r2_3) * x + r2_2) * x + r2_1) * x + r2_0) /+ ((((x + r2_8) * x + r2_7) * x + r2_6) * x + r2_5)+ | x < 12 = ((((r3_4 * x + r3_3) * x + r3_2) * x + r3_1) * x + r3_0) /+ ((((x + r3_8) * x + r3_7) * x + r3_6) * x + r3_5)+ | x > 3e6 = k+ | otherwise = k + x1 *+ ((r4_2 * x2 + r4_1) * x2 + r4_0) /+ ((x2 + r4_4) * x2 + r4_3)+ where+ (a , b , c)+ | x < 0.5 = (-y , x + 1 , x)+ | otherwise = (0 , x , x - 1)++ y = log x+ k = x * (y-1) - 0.5 * y + alr2pi+ alr2pi = 0.918938533204673++ x1 = 1 / x+ x2 = x1 * x1++ r1_0 = -2.66685511495; r1_1 = -24.4387534237; r1_2 = -21.9698958928+ r1_3 = 11.1667541262; r1_4 = 3.13060547623; r1_5 = 0.607771387771+ r1_6 = 11.9400905721; r1_7 = 31.4690115749; r1_8 = 15.2346874070++ r2_0 = -78.3359299449; r2_1 = -142.046296688; r2_2 = 137.519416416+ r2_3 = 78.6994924154; r2_4 = 4.16438922228; r2_5 = 47.0668766060+ r2_6 = 313.399215894; r2_7 = 263.505074721; r2_8 = 43.3400022514++ r3_0 = -2.12159572323e5; r3_1 = 2.30661510616e5; r3_2 = 2.74647644705e4+ r3_3 = -4.02621119975e4; r3_4 = -2.29660729780e3; r3_5 = -1.16328495004e5+ r3_6 = -1.46025937511e5; r3_7 = -2.42357409629e4; r3_8 = -5.70691009324e2++ r4_0 = 0.279195317918525; r4_1 = 0.4917317610505968;+ r4_2 = 0.0692910599291889; r4_3 = 3.350343815022304+ r4_4 = 6.012459259764103
Numeric/SpecFunctions/Internal.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE BangPatterns, ScopedTypeVariables #-}+{-# LANGUAGE BangPatterns, ScopedTypeVariables, ForeignFunctionInterface #-} -- | -- Module : Numeric.SpecFunctions.Internal -- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan@@ -11,16 +11,20 @@ -- Internal module with implementation of special functions. module Numeric.SpecFunctions.Internal where +import Control.Applicative import Data.Bits ((.&.), (.|.), shiftR) import Data.Int (Int64)-import qualified Data.Number.Erf as Erf (erfc,erf)+import Data.Word (Word) import qualified Data.Vector.Unboxed as U+import Data.Vector.Unboxed ((!)) import Numeric.Polynomial.Chebyshev (chebyshevBroucke)-import Numeric.Polynomial (evaluateEvenPolynomialL,evaluateOddPolynomialL)+import Numeric.Polynomial (evaluatePolynomialL,evaluateEvenPolynomialL,evaluateOddPolynomialL)+import Numeric.RootFinding (Root(..), newtonRaphson)+import Numeric.Series (sumPowerSeries,enumSequenceFrom,scanSequence,evalContFractionB) import Numeric.MathFunctions.Constants ( m_epsilon, m_NaN, m_neg_inf, m_pos_inf- , m_sqrt_2_pi, m_ln_sqrt_2_pi, m_sqrt_2- , m_eulerMascheroni+ , m_sqrt_2_pi, m_ln_sqrt_2_pi, m_eulerMascheroni+ , m_min_log, m_tiny ) import Text.Printf @@ -31,23 +35,46 @@ -- | Error function. ----- > erf -∞ = -1--- > erf 0 = 0--- > erf +∞ = 1+-- \[+-- \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} \exp(-t^2) dt+-- \]+--+-- Function limits are:+--+-- \[+-- \begin{aligned}+-- &\operatorname{erf}(-\infty) &=& -1 \\+-- &\operatorname{erf}(0) &=& \phantom{-}\,0 \\+-- &\operatorname{erf}(+\infty) &=& \phantom{-}\,1 \\+-- \end{aligned}+-- \] erf :: Double -> Double {-# INLINE erf #-}-erf = Erf.erf+erf = c_erf -- | Complementary error function. ----- > erfc -∞ = 2--- > erfc 0 = 1--- > errc +∞ = 0+-- \[+-- \operatorname{erfc}(x) = 1 - \operatorname{erf}(x)+-- \]+--+-- Function limits are:+--+-- \[+-- \begin{aligned}+-- &\operatorname{erf}(-\infty) &=&\, 2 \\+-- &\operatorname{erf}(0) &=&\, 1 \\+-- &\operatorname{erf}(+\infty) &=&\, 0 \\+-- \end{aligned}+-- \] erfc :: Double -> Double {-# INLINE erfc #-}-erfc = Erf.erfc+erfc = c_erfc +foreign import ccall "erf" c_erf :: Double -> Double+foreign import ccall "erfc" c_erfc :: Double -> Double + -- | Inverse of 'erf'. invErf :: Double -- ^ /p/ ∈ [-1,1] -> Double@@ -81,103 +108,59 @@ -- Gamma function ---------------------------------------------------------------- --- Adapted from http://people.sc.fsu.edu/~burkardt/f_src/asa245/asa245.html+data L = L {-# UNPACK #-} !Double {-# UNPACK #-} !Double --- | Compute the logarithm of the gamma function Γ(/x/). Uses--- Algorithm AS 245 by Macleod.+-- | Compute the logarithm of the gamma function, Γ(/x/). ----- Gives an accuracy of 10-12 significant decimal digits, except--- for small regions around /x/ = 1 and /x/ = 2, where the function--- goes to zero. For greater accuracy, use 'logGammaL'.+-- \[+-- \Gamma(x) = \int_0^{\infty}t^{x-1}e^{-t}\,dt = (x - 1)!+-- \] ----- Returns ∞ if the input is outside of the range (0 < /x/ ≤ 1e305).+-- This implementation uses Lanczos approximation. It gives 14 or more+-- significant decimal digits, except around /x/ = 1 and /x/ = 2,+-- where the function goes to zero.+--+-- Returns ∞ if the input is outside of the range (0 < /x/+-- ≤ 1e305). logGamma :: Double -> Double logGamma x- | x <= 0 = m_pos_inf- -- Handle positive infinity. logGamma overflows before 1e308 so- -- it's safe- | x > 1e308 = m_pos_inf- -- Normal cases- | x < 1.5 = a + c *- ((((r1_4 * b + r1_3) * b + r1_2) * b + r1_1) * b + r1_0) /- ((((b + r1_8) * b + r1_7) * b + r1_6) * b + r1_5)- | x < 4 = (x - 2) *- ((((r2_4 * x + r2_3) * x + r2_2) * x + r2_1) * x + r2_0) /- ((((x + r2_8) * x + r2_7) * x + r2_6) * x + r2_5)- | x < 12 = ((((r3_4 * x + r3_3) * x + r3_2) * x + r3_1) * x + r3_0) /- ((((x + r3_8) * x + r3_7) * x + r3_6) * x + r3_5)- | x > 3e6 = k- | otherwise = k + x1 *- ((r4_2 * x2 + r4_1) * x2 + r4_0) /- ((x2 + r4_4) * x2 + r4_3)+ | x <= 0 = m_pos_inf+ | x < 1 = lanczos (1+x) - log x+ | otherwise = lanczos x where- (a , b , c)- | x < 0.5 = (-y , x + 1 , x)- | otherwise = (0 , x , x - 1)-- y = log x- k = x * (y-1) - 0.5 * y + alr2pi- alr2pi = 0.918938533204673-- x1 = 1 / x- x2 = x1 * x1-- r1_0 = -2.66685511495; r1_1 = -24.4387534237; r1_2 = -21.9698958928- r1_3 = 11.1667541262; r1_4 = 3.13060547623; r1_5 = 0.607771387771- r1_6 = 11.9400905721; r1_7 = 31.4690115749; r1_8 = 15.2346874070-- r2_0 = -78.3359299449; r2_1 = -142.046296688; r2_2 = 137.519416416- r2_3 = 78.6994924154; r2_4 = 4.16438922228; r2_5 = 47.0668766060- r2_6 = 313.399215894; r2_7 = 263.505074721; r2_8 = 43.3400022514-- r3_0 = -2.12159572323e5; r3_1 = 2.30661510616e5; r3_2 = 2.74647644705e4- r3_3 = -4.02621119975e4; r3_4 = -2.29660729780e3; r3_5 = -1.16328495004e5- r3_6 = -1.46025937511e5; r3_7 = -2.42357409629e4; r3_8 = -5.70691009324e2-- r4_0 = 0.279195317918525; r4_1 = 0.4917317610505968;- r4_2 = 0.0692910599291889; r4_3 = 3.350343815022304- r4_4 = 6.012459259764103----data L = L {-# UNPACK #-} !Double {-# UNPACK #-} !Double+ -- Evaluate Lanczos approximation for γ=6+ lanczos z = fini+ $ U.foldl' go (L 0 (z+7)) a+ where+ fini (L l _) = log (l+a0) + log m_sqrt_2_pi - z65 + (z-0.5) * log z65+ go (L l t) k = L (l + k / t) (t-1)+ z65 = z + 6.5+ -- Coefficients for Lanczos approximation+ a0 = 0.9999999999995183+ a = U.fromList [ 0.1659470187408462e-06+ , 0.9934937113930748e-05+ , -0.1385710331296526+ , 12.50734324009056+ , -176.6150291498386+ , 771.3234287757674+ , -1259.139216722289+ , 676.5203681218835+ ] --- | Compute the logarithm of the gamma function, Γ(/x/). Uses a--- Lanczos approximation.------ This function is slower than 'logGamma', but gives 14 or more--- significant decimal digits of accuracy, except around /x/ = 1 and--- /x/ = 2, where the function goes to zero.------ Returns ∞ if the input is outside of the range (0 < /x/--- ≤ 1e305).+-- | Synonym for 'logGamma'. Retained for compatibility logGammaL :: Double -> Double-logGammaL x- | x <= 0 = m_pos_inf- -- Lanroz approximation loses precision for small arguments- | x <= 1e-3 = logGamma x- | otherwise = fini . U.foldl' go (L 0 (x+7)) $ a- where fini (L l _) = log (l+a0) + log m_sqrt_2_pi - x65 + (x-0.5) * log x65- go (L l t) k = L (l + k / t) (t-1)- x65 = x + 6.5- a0 = 0.9999999999995183- a = U.fromList [ 0.1659470187408462e-06- , 0.9934937113930748e-05- , -0.1385710331296526- , 12.50734324009056- , -176.6150291498386- , 771.3234287757674- , -1259.139216722289- , 676.5203681218835- ]-+logGammaL = logGamma --- | Compute the log gamma correction factor for @x@ ≥ 10. This--- correction factor is suitable for an alternate (but less--- numerically accurate) definition of 'logGamma':+-- |+-- Compute the log gamma correction factor for Stirling+-- approximation for @x@ ≥ 10. This correction factor is+-- suitable for an alternate (but less numerically accurate)+-- definition of 'logGamma': ----- >lgg x = 0.5 * log(2*pi) + (x-0.5) * log x - x + logGammaCorrection x+-- \[+-- \log\Gamma(x) = \frac{1}{2}\log(2\pi) + (x-\frac{1}{2})\log x - x + \operatorname{logGammaCorrection}(x)+-- \] logGammaCorrection :: Double -> Double logGammaCorrection x | x < 10 = m_NaN@@ -199,83 +182,109 @@ -- | Compute the normalized lower incomplete gamma function--- γ(/s/,/x/). Normalization means that--- γ(/s/,∞)=1. Uses Algorithm AS 239 by Shea.-incompleteGamma :: Double -- ^ /s/ ∈ (0,∞)+-- γ(/z/,/x/). Normalization means that γ(/z/,∞)=1+--+-- \[+-- \gamma(z,x) = \frac{1}{\Gamma(z)}\int_0^{x}t^{z-1}e^{-t}\,dt+-- \]+--+-- Uses Algorithm AS 239 by Shea.+incompleteGamma :: Double -- ^ /z/ ∈ (0,∞) -> Double -- ^ /x/ ∈ (0,∞) -> Double-incompleteGamma p x- | isNaN p || isNaN x = m_NaN- | x < 0 || p <= 0 = m_pos_inf- | x == 0 = 0- -- For very large `p' normal approximation gives <1e-10 error- | p >= 2e5 = norm (3 * sqrt p * ((x/p) ** (1/3) + 1/(9*p) - 1))- | p >= 500 = approx- -- Dubious approximation- | x >= 1e8 = 1- | x <= 1 || x < p = let a = p * log x - x - logGamma (p + 1)- g = a + log (pearson p 1 1)- in if g > limit then exp g else 0- | otherwise = let g = p * log x - x - logGamma p + log cf- in if g > limit then 1 - exp g else 1+-- Notation used:+-- + P(a,x) - regularized lower incomplete gamma+-- + Q(a,x) - regularized upper incomplete gamma+incompleteGamma a x+ | a <= 0 || x < 0 = error+ $ "incompleteGamma: Domain error z=" ++ show a ++ " x=" ++ show x+ | x == 0 = 0+ | x == m_pos_inf = 1+ -- For very small x we use following expansion for P:+ --+ -- See http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/+ | x < sqrt m_epsilon && a > 1+ = x**a / a / exp (logGammaL a) * (1 - a*x / (a + 1))+ | x < 0.5 = case () of+ _| (-0.4)/log x < a -> taylorSeriesP+ | otherwise -> taylorSeriesComplQ+ | x < 1.1 = case () of+ _| 0.75*x < a -> taylorSeriesP+ | otherwise -> taylorSeriesComplQ+ | a > 20 && useTemme = uniformExpansion+ | x - (1 / (3 * x)) < a = taylorSeriesP+ | otherwise = contFraction where- -- CDF for standard normal distributions- norm a = 0.5 * erfc (- a / m_sqrt_2)- -- For large values of `p' we use 18-point Gauss-Legendre- -- integration.- approx- | ans > 0 = 1 - ans- | otherwise = -ans- where- -- Set upper limit for integration- xu | x > p1 = (p1 + 11.5*sqrtP1) `max` (x + 6*sqrtP1)- | otherwise = max 0 $ (p1 - 7.5*sqrtP1) `min` (x - 5*sqrtP1)- s = U.sum $ U.zipWith go coefY coefW- go y w = let t = x + (xu - x)*y- in w * exp( -(t-p1) + p1*(log t - lnP1) )- ans = s * (xu - x) * exp( p1 * (lnP1 - 1) - logGamma p)- --- p1 = p - 1- lnP1 = log p1- sqrtP1 = sqrt p1+ mu = (x - a) / a+ useTemme = (a > 200 && 20/a > mu*mu)+ || (abs mu < 0.4)+ -- Gautschi's algorithm. --- pearson !a !c !g- | c' <= tolerance = g'- | otherwise = pearson a' c' g'- where a' = a + 1- c' = c * x / a'- g' = g + c'- cf = let a = 1 - p- b = a + x + 1- p3 = x + 1- p4 = x * b- in contFrac a b 0 1 x p3 p4 (p3/p4)- contFrac !a !b !c !p1 !p2 !p3 !p4 !g- | abs (g - rn) <= min tolerance (tolerance * rn) = g- | otherwise = contFrac a' b' c' (f p3) (f p4) (f p5) (f p6) rn- where a' = a + 1- b' = b + 2- c' = c + 1- an = a' * c'- p5 = b' * p3 - an * p1- p6 = b' * p4 - an * p2- rn = p5 / p6- f n | abs p5 > overflow = n / overflow- | otherwise = n- limit = -88- tolerance = 1e-14- overflow = 1e37+ -- Evaluate series for P(a,x). See [Temme1994] Eq. 5.5+ --+ -- FIXME: Term `exp (log x * z - x - logGamma (z+1))` doesn't give full precision+ taylorSeriesP+ = sumPowerSeries x (scanSequence (/) 1 $ enumSequenceFrom (a+1))+ * exp (log x * a - x - logGamma (a+1))+ -- Series for 1-Q(a,x). See [Temme1994] Eq. 5.5+ taylorSeriesComplQ+ = sumPowerSeries (-x) (scanSequence (/) 1 (enumSequenceFrom 1) / enumSequenceFrom a)+ * x**a / exp(logGammaL a)+ -- Legendre continued fractions+ contFraction = 1 - ( exp ( log x * a - x - logGamma a )+ / evalContFractionB frac+ )+ where+ frac = (\k -> (k*(a-k), x - a + 2*k + 1)) <$> enumSequenceFrom 0+ -- Evaluation based on uniform expansions. See [Temme1994] 5.2+ uniformExpansion =+ let -- Coefficients f_m in paper+ fm :: U.Vector Double+ fm = U.fromList [ 1.00000000000000000000e+00+ ,-3.33333333333333370341e-01+ , 8.33333333333333287074e-02+ ,-1.48148148148148153802e-02+ , 1.15740740740740734316e-03+ , 3.52733686067019369930e-04+ ,-1.78755144032921825352e-04+ , 3.91926317852243766954e-05+ ,-2.18544851067999240532e-06+ ,-1.85406221071515996597e-06+ , 8.29671134095308545622e-07+ ,-1.76659527368260808474e-07+ , 6.70785354340149841119e-09+ , 1.02618097842403069078e-08+ ,-4.38203601845335376897e-09+ , 9.14769958223679020897e-10+ ,-2.55141939949462514346e-11+ ,-5.83077213255042560744e-11+ , 2.43619480206674150369e-11+ ,-5.02766928011417632057e-12+ , 1.10043920319561347525e-13+ , 3.37176326240098513631e-13+ ]+ y = - log1pmx mu+ eta = sqrt (2 * y) * signum mu+ -- Evaluate S_α (Eq. 5.9)+ loop !_ !_ u 0 = u+ loop bm1 bm0 u i = let t = (fm ! i) + (fromIntegral i + 1)*bm1 / a+ u' = eta * u + t+ in loop bm0 t u' (i-1)+ s_a = let n = U.length fm+ in loop (fm ! (n-1)) (fm ! (n-2)) 0 (n-3)+ / exp (logGammaCorrection a)+ in 1/2 * erfc(-eta*sqrt(a/2)) - exp(-(a*y)) / sqrt (2*pi*a) * s_a -- Adapted from Numerical Recipes §6.2.1 -- | Inverse incomplete gamma function. It's approximately inverse of--- 'incompleteGamma' for the same /s/. So following equality+-- 'incompleteGamma' for the same /z/. So following equality -- approximately holds: ----- > invIncompleteGamma s . incompleteGamma s = id-invIncompleteGamma :: Double -- ^ /s/ ∈ (0,∞)+-- > invIncompleteGamma z . incompleteGamma z ≈ id+invIncompleteGamma :: Double -- ^ /z/ ∈ (0,∞) -> Double -- ^ /p/ ∈ [0,1] -> Double invIncompleteGamma a p@@ -341,13 +350,27 @@ ---------------------------------------------------------------- -- | Compute the natural logarithm of the beta function.-logBeta :: Double -> Double -> Double+--+-- \[+-- B(a,b) = \int_0^1 t^{a-1}(1-t)^{1-b}\,dt = \frac{\Gamma{a}\Gamma{b}}{\Gamma{a+b}}+-- \]+logBeta+ :: Double -- ^ /a/ > 0+ -> Double -- ^ /b/ > 0+ -> Double logBeta a b | p < 0 = m_NaN | p == 0 = m_pos_inf- | p >= 10 = log q * (-0.5) + m_ln_sqrt_2_pi + logGammaCorrection p + c +- (p - 0.5) * log ppq + q * log1p(-ppq)- | q >= 10 = logGamma p + c + p - p * log pq + (q - 0.5) * log1p(-ppq)+ | p >= 10 = log q * (-0.5)+ + m_ln_sqrt_2_pi+ + logGammaCorrection p+ + c+ + (p - 0.5) * log ppq + q * log1p(-ppq)+ | q >= 10 = logGamma p+ + c+ + p+ - p * log pq+ + (q - 0.5) * log1p(-ppq) | otherwise = logGamma p + logGamma q - logGamma pq where p = min a b@@ -356,11 +379,16 @@ pq = p + q c = logGammaCorrection q - logGammaCorrection pq --- | Regularized incomplete beta function. Uses algorithm AS63 by--- Majumder and Bhattachrjee and quadrature approximation for large--- /p/ and /q/.-incompleteBeta :: Double -- ^ /p/ > 0- -> Double -- ^ /q/ > 0+-- | Regularized incomplete beta function.+--+-- \[+-- I(x;a,b) = \frac{1}{B(a,b)} \int_0^x t^{a-1}(1-t)^{1-b}\,dt+-- \]+--+-- Uses algorithm AS63 by Majumder and Bhattachrjee and quadrature+-- approximation for large /p/ and /q/.+incompleteBeta :: Double -- ^ /a/ > 0+ -> Double -- ^ /b/ > 0 -> Double -- ^ /x/, must lie in [0,1] range -> Double incompleteBeta p q = incompleteBeta_ (logBeta p q) p q@@ -368,15 +396,15 @@ -- | Regularized incomplete beta function. Same as 'incompleteBeta' -- but also takes logarithm of beta function as parameter. incompleteBeta_ :: Double -- ^ logarithm of beta function for given /p/ and /q/- -> Double -- ^ /p/ > 0- -> Double -- ^ /q/ > 0+ -> Double -- ^ /a/ > 0+ -> Double -- ^ /b/ > 0 -> Double -- ^ /x/, must lie in [0,1] range -> Double incompleteBeta_ beta p q x | p <= 0 || q <= 0 = modErr $ printf "incompleteBeta_: p <= 0 || q <= 0. p=%g q=%g x=%g" p q x | x < 0 || x > 1 || isNaN x =- modErr $ printf "incompletBeta_: x out of [0,1] range. p=%g q=%g x=%g" p q x+ modErr $ printf "incompleteBeta_: x out of [0,1] range. p=%g q=%g x=%g" p q x | x == 0 || x == 1 = x | p >= (p+q) * x = incompleteBetaWorker beta p q x | otherwise = 1 - incompleteBetaWorker beta q p (1 - x)@@ -394,8 +422,8 @@ p1 = p - 1 q1 = q - 1 mu = p / (p + q)- lnmu = log mu- lnmuc = log (1 - mu)+ lnmu = log mu+ lnmuc = log1p (-mu) -- Upper limit for integration xu = max 0 $ min (mu - 10*t) (x - 5*t) where@@ -419,9 +447,22 @@ -- Constants eps = 1e-15 cx = 1 - x- -- Loop+ -- Common multiplies for expansion. Accurate calculation is a bit+ -- tricky. Performing calculation in log-domain leads to slight+ -- loss of precision for small x, while using ** prone to+ -- underflows.+ --+ -- If either beta function of x**p·(1-x)**(q-1) underflows we+ -- switch to log domain. It could waste work but there's no easy+ -- switch criterion.+ factor+ | beta < m_min_log || prod < m_tiny = exp( p * log x + (q - 1) * log cx - beta)+ | otherwise = prod / exp beta+ where+ prod = x**p * cx**(q - 1)+ -- Soper's expansion of incomplete beta function loop !psq (ns :: Int) ai term betain- | done = betain' * exp( p * log x + (q - 1) * log cx - beta) / p+ | done = betain' * factor / p | otherwise = loop psq' (ns - 1) (ai + 1) term' betain' where -- New values@@ -439,9 +480,9 @@ -- | Compute inverse of regularized incomplete beta function. Uses -- initial approximation from AS109, AS64 and Halley method to solve -- equation.-invIncompleteBeta :: Double -- ^ /p/ > 0- -> Double -- ^ /q/ > 0- -> Double -- ^ /a/ ∈ [0,1]+invIncompleteBeta :: Double -- ^ /a/ > 0+ -> Double -- ^ /b/ > 0+ -> Double -- ^ /x/ ∈ [0,1] -> Double invIncompleteBeta p q a | p <= 0 || q <= 0 =@@ -455,7 +496,7 @@ invIncompleteBetaWorker :: Double -> Double -> Double -> Double -> Double -- NOTE: p <= 0.5.-invIncompleteBetaWorker beta a b p = loop (0::Int) guess+invIncompleteBetaWorker beta a b p = loop (0::Int) (invIncBetaGuess beta a b p) where a1 = a - 1 b1 = b - 1@@ -478,64 +519,180 @@ where -- Calculate Halley step. f = incompleteBeta_ beta a b x - p- f' = exp $ a1 * log x + b1 * log (1 - x) - beta+ f' = exp $ a1 * log x + b1 * log1p (-x) - beta u = f / f'- dx = u / (1 - 0.5 * min 1 (u * (a1 / x - b1 / (1 - x))))+ -- We bound Halley correction to Newton-Raphson to (-1,1) range+ corr | d > 1 = 1+ | d < -1 = -1+ | isNaN d = 0+ | otherwise = d+ where+ d = u * (a1 / x - b1 / (1 - x))+ dx = u / (1 - 0.5 * corr) -- Next approximation. If Halley step leads us out of [0,1] -- range we revert to bisection. x' | z < 0 = x / 2 | z > 1 = (x + 1) / 2 | otherwise = z where z = x - dx- -- Calculate initial guess. Approximations from AS64, AS109 and- -- Numerical recipes are used.- --- -- Equations are referred to by name of paper and number e.g. [AS64 2]- -- In AS64 papers equations are not numbered so they are refered- -- to by number of appearance starting from definition of- -- incomplete beta.- guess- -- In this region we use approximation from AS109 (Carter- -- approximation). It's reasonably good (2 iterations on- -- average)- | a > 1 && b > 1 =- let r = (y*y - 3) / 6- s = 1 / (2*a - 1)- t = 1 / (2*b - 1)- h = 2 / (s + t)- w = y * sqrt(h + r) / h - (t - s) * (r + 5/6 - 2 / (3 * h))- in a / (a + b * exp(2 * w))- -- Otherwise we revert to approximation from AS64 derived from- -- [AS64 2] when it's applicable.- --- -- It slightly reduces average number of iterations when `a' and- -- `b' have different magnitudes.- | chi2 > 0 && ratio > 1 = 1 - 2 / (ratio + 1)- -- If all else fails we use approximation from "Numerical- -- Recipes". It's very similar to approximations [AS64 4,5] but- -- it never goes out of [0,1] interval.- | otherwise = case () of- _| p < t / w -> (a * p * w) ** (1/a)- | otherwise -> 1 - (b * (1 - p) * w) ** (1/b)- where- lna = log $ a / (a+b)- lnb = log $ b / (a+b)- t = exp( a * lna ) / a- u = exp( b * lnb ) / b- w = t + u+++-- Calculate initial guess for inverse incomplete beta function.+invIncBetaGuess :: Double -> Double -> Double -> Double -> Double+-- Calculate initial guess. for solving equation for inverse incomplete beta.+-- It's really hodgepodge of different approximations accumulated over years.+--+-- Equations are referred to by name of paper and number e.g. [AS64 2]+-- In AS64 papers equations are not numbered so they are refered to by+-- number of appearance starting from definition of incomplete beta.+invIncBetaGuess beta a b p+ -- If both a and b are less than 1 incomplete beta have inflection+ -- point.+ --+ -- > x = (1 - a) / (2 - a - b)+ --+ -- We approximate incomplete beta by neglecting one of factors under+ -- integral and then rescaling result of integration into [0,1]+ -- range.+ | a < 1 && b < 1 =+ let x_infl = (1 - a) / (2 - a - b)+ p_infl = incompleteBeta a b x_infl+ x | p < p_infl = let xg = (a * p * exp beta) ** (1/a) in xg / (1+xg)+ | otherwise = let xg = (b * (1-p) * exp beta) ** (1/b) in 1 - xg/(1+xg)+ in x+ -- If both a and b larger or equal that 1 but not too big we use+ -- same approximation as above but calculate it a bit differently+ | a+b <= 6 && a>=1 && b>=1 =+ let x_infl = (a - 1) / (a + b - 2)+ p_infl = incompleteBeta a b x_infl+ x | p < p_infl = exp ((log(p * a) + beta) / a)+ | otherwise = 1 - exp((log((1-p) * b) + beta) / b)+ in x+ -- For small a and not too big b we use approximation from boost.+ | b < 5 && a < 1 =+ let x | p**(1/a) < 0.5 = (p * a * exp beta) ** (1/a)+ | otherwise = 1 - (1 - p ** (b * exp beta))**(1/b)+ in x+ -- When a>>b and both are large approximation from [Temme1992],+ -- section 4 "the incomplete gamma function case" used. In this+ -- region it greatly improves over other approximation (AS109, AS64,+ -- "Numerical Recipes")+ --+ -- FIXME: It could be used when b>>a too but it require inverse of+ -- upper incomplete gamma to be precise enough. In current+ -- implementation it loses precision in horrible way (40+ -- order of magnitude off for sufficiently small p)+ | a+b > 5 && a/b > 4 =+ let -- Calculate initial approximation to eta using eq 4.1+ eta0 = invIncompleteGamma b (1-p) / a+ mu = b / a -- Eq. 4.3+ -- A lot of helpers for calculation of+ w = sqrt(1 + mu) -- Eq. 4.9+ w_2 = w * w+ w_3 = w_2 * w+ w_4 = w_2 * w_2+ w_5 = w_3 * w_2+ w_6 = w_3 * w_3+ w_7 = w_4 * w_3+ w_8 = w_4 * w_4+ w_9 = w_5 * w_4+ w_10 = w_5 * w_5+ d = eta0 - mu+ d_2 = d * d+ d_3 = d_2 * d+ d_4 = d_2 * d_2+ w1 = w + 1+ w1_2 = w1 * w1+ w1_3 = w1 * w1_2+ w1_4 = w1_2 * w1_2+ -- Evaluation of eq 4.10+ e1 = (w + 2) * (w - 1) / (3 * w)+ + (w_3 + 9 * w_2 + 21 * w + 5) * d+ / (36 * w_2 * w1)+ - (w_4 - 13 * w_3 + 69 * w_2 + 167 * w + 46) * d_2+ / (1620 * w1_2 * w_3)+ - (7 * w_5 + 21 * w_4 + 70 * w_3 + 26 * w_2 - 93 * w - 31) * d_3+ / (6480 * w1_3 * w_4)+ - (75 * w_6 + 202 * w_5 + 188 * w_4 - 888 * w_3 - 1345 * w_2 + 118 * w + 138) * d_4+ / (272160 * w1_4 * w_5)+ e2 = (28 * w_4 + 131 * w_3 + 402 * w_2 + 581 * w + 208) * (w - 1)+ / (1620 * w1 * w_3)+ - (35 * w_6 - 154 * w_5 - 623 * w_4 - 1636 * w_3 - 3983 * w_2 - 3514 * w - 925) * d+ / (12960 * w1_2 * w_4)+ - ( 2132 * w_7 + 7915 * w_6 + 16821 * w_5 + 35066 * w_4 + 87490 * w_3+ + 141183 * w_2 + 95993 * w + 21640+ ) * d_2+ / (816480 * w_5 * w1_3)+ - ( 11053 * w_8 + 53308 * w_7 + 117010 * w_6 + 163924 * w_5 + 116188 * w_4+ - 258428 * w_3 - 677042 * w_2 - 481940 * w - 105497+ ) * d_3+ / (14696640 * w1_4 * w_6)+ e3 = -( (3592 * w_7 + 8375 * w_6 - 1323 * w_5 - 29198 * w_4 - 89578 * w_3+ - 154413 * w_2 - 116063 * w - 29632+ ) * (w - 1)+ )+ / (816480 * w_5 * w1_2)+ - ( 442043 * w_9 + 2054169 * w_8 + 3803094 * w_7 + 3470754 * w_6 + 2141568 * w_5+ - 2393568 * w_4 - 19904934 * w_3 - 34714674 * w_2 - 23128299 * w - 5253353+ ) * d+ / (146966400 * w_6 * w1_3)+ - ( 116932 * w_10 + 819281 * w_9 + 2378172 * w_8 + 4341330 * w_7 + 6806004 * w_6+ + 10622748 * w_5 + 18739500 * w_4 + 30651894 * w_3 + 30869976 * w_2+ + 15431867 * w + 2919016+ ) * d_2+ / (146966400 * w1_4 * w_7)+ eta = evaluatePolynomialL (1/a) [eta0, e1, e2, e3]+ -- Now we solve eq 4.2 to recover x using Newton iterations+ u = eta - mu * log eta + (1 + mu) * log(1 + mu) - mu+ cross = 1 / (1 + mu);+ lower = if eta < mu then cross else 0+ upper = if eta < mu then 1 else cross+ x_guess = (lower + upper) / 2+ func x = ( u + log x + mu*log(1 - x)+ , 1/x - mu/(1-x)+ )+ Root x0 = newtonRaphson 1e-8 (lower, x_guess, upper) func+ in x0+ -- For large a and b approximation from AS109 (Carter+ -- approximation). It's reasonably good in this region+ | a > 1 && b > 1 =+ let r = (y*y - 3) / 6+ s = 1 / (2*a - 1)+ t = 1 / (2*b - 1)+ h = 2 / (s + t)+ w = y * sqrt(h + r) / h - (t - s) * (r + 5/6 - 2 / (3 * h))+ in a / (a + b * exp(2 * w))+ -- Otherwise we revert to approximation from AS64 derived from+ -- [AS64 2] when it's applicable.+ --+ -- It slightly reduces average number of iterations when `a' and+ -- `b' have different magnitudes.+ | chi2 > 0 && ratio > 1 = 1 - 2 / (ratio + 1)+ -- If all else fails we use approximation from "Numerical+ -- Recipes". It's very similar to approximations [AS64 4,5] but+ -- it never goes out of [0,1] interval.+ | otherwise = case () of+ _| p < t / w -> (a * p * w) ** (1/a)+ | otherwise -> 1 - (b * (1 - p) * w) ** (1/b)+ where+ lna = log $ a / (a+b)+ lnb = log $ b / (a+b)+ t = exp( a * lna ) / a+ u = exp( b * lnb ) / b+ w = t + u+ where+ -- Formula [AS64 2]+ ratio = (4*a + 2*b - 2) / chi2+ -- Quantile of chi-squared distribution. Formula [AS64 3].+ chi2 = 2 * b * (1 - t + y * sqrt t) ** 3 where- -- Formula [2]- ratio = (4*a + 2*b - 2) / chi2- -- Quantile of chi-squared distribution. Formula [3].- chi2 = 2 * b * (1 - t + y * sqrt t) ** 3- where- t = 1 / (9 * b)- -- `y' is Hasting's approximation of p'th quantile of standard- -- normal distribution.- y = r - ( 2.30753 + 0.27061 * r )- / ( 1.0 + ( 0.99229 + 0.04481 * r ) * r )- where- r = sqrt $ - 2 * log p+ t = 1 / (9 * b)+ -- `y' is Hasting's approximation of p'th quantile of standard+ -- normal distribution.+ y = r - ( 2.30753 + 0.27061 * r )+ / ( 1.0 + ( 0.99229 + 0.04481 * r ) * r )+ where+ r = sqrt $ - 2 * log p @@ -603,7 +760,25 @@ -0.10324619158271569595141333961932e-15 ] +-- | Compute log(1+x)-x:+log1pmx :: Double -> Double+log1pmx x+ | x < -1 = error "Domain error"+ | x == -1 = m_neg_inf+ | ax > 0.95 = log(1 + x) - x+ | ax < m_epsilon = -(x * x) /2+ | otherwise = - x * x * sumPowerSeries (-x) (recip <$> enumSequenceFrom 2)+ where+ ax = abs x +-- | Compute @exp x - 1@ without loss of accuracy for x near zero.+expm1 :: Double -> Double+expm1 = c_expm1++foreign import ccall "expm1" c_expm1 :: Double -> Double+++ -- | /O(log n)/ Compute the logarithm in base 2 of the given value. log2 :: Int -> Int log2 v0@@ -615,7 +790,9 @@ in go (i-1) (r .|. si) (v `shiftR` si) | otherwise = go (i-1) r v b = U.unsafeIndex bv- !bv = U.fromList [0x2, 0xc, 0xf0, 0xff00, 0xffff0000, 0xffffffff00000000]+ !bv = U.fromList [ 0x02, 0x0c, 0xf0, 0xff00+ , fromIntegral (0xffff0000 :: Word)+ , fromIntegral (0xffffffff00000000 :: Word)] !sv = U.fromList [1,2,4,8,16,32] @@ -742,9 +919,14 @@ max64 = fromIntegral (maxBound :: Int64) round64 x = round x :: Int64 --- | Compute ψ0(/x/), the first logarithmic derivative of the gamma--- function. Uses Algorithm AS 103 by Bernardo, based on Minka's C--- implementation.+-- | Compute ψ(/x/), the first logarithmic derivative of the gamma+-- function.+--+-- \[+-- \psi(x) = \frac{d}{dx} \ln \left(\Gamma(x)\right) = \frac{\Gamma'(x)}{\Gamma(x)}+-- \]+--+-- Uses Algorithm AS 103 by Bernardo, based on Minka's C implementation. digamma :: Double -> Double digamma x | isNaN x || isInfinite x = m_NaN
Numeric/Sum.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE BangPatterns, CPP, DeriveDataTypeable, FlexibleContexts,+{-# LANGUAGE BangPatterns, DeriveDataTypeable, FlexibleContexts, MultiParamTypeClasses, TemplateHaskell, TypeFamilies #-} {-# OPTIONS_GHC -fno-warn-name-shadowing #-} -- |
changelog.md view
@@ -1,3 +1,36 @@+Changes in 0.2.0.0++ * `logGamma` now uses Lancsoz approximation and same as `logGammaL`. Old+ implementation of `logGamma` moved to `Numeric.SpecFunctions.Extra.logGammaAS245`.++ * Precision of `logGamma` for z<1 improved.++ * New much more precise implementation for `incompleteGamma`++ * Dependency on `erf` pacakge dropped. `erf` and `erfc` just do direct calls+ to C.++ * `Numeric.SpecFunctions.expm1` added++ * `Numeric.SpecFunctions.log1pmx` added.++ * `logGammaCorrection` exported in `Numeric.SpecFunctions.Extra`.++ * Module `Numeric.Series` added for working with infinite sequences, series+ summation and evaluation of continued fractions.++ * Module `statistics: Statistics.Math.RootFinding` copied to+ `Numeric.RootFinding`. Instances for `binary` and `aeson` dropped.++ * Root-finding using Newton-Raphson added++ * `Numeric.MathFunctions.Comparison.ulpDelta` added. It calculates signed+ distance between two doubles.++ * Other bug fixes.+++ Changes in 0.1.7.0 * Module `statistics: Statistics.Function.Comparison` moved to
math-functions.cabal view
@@ -1,6 +1,6 @@ name: math-functions-version: 0.1.7.0-cabal-version: >= 1.8+version: 0.2.0.0+cabal-version: >= 1.10 license: BSD3 license-file: LICENSE author: Bryan O'Sullivan <bos@serpentine.com>,@@ -26,10 +26,20 @@ doc/sinc.hs library- ghc-options: -Wall+ default-language: Haskell2010+ other-extensions:+ BangPatterns+ CPP+ DeriveDataTypeable+ FlexibleContexts+ MultiParamTypeClasses+ ScopedTypeVariables+ TemplateHaskell+ TypeFamilies++ ghc-options: -Wall -O2 build-depends: base >=3 && <5, deepseq,- erf >= 2, vector >= 0.7, primitive, vector-th-unbox@@ -38,13 +48,18 @@ Numeric.MathFunctions.Comparison Numeric.Polynomial Numeric.Polynomial.Chebyshev+ Numeric.RootFinding Numeric.SpecFunctions Numeric.SpecFunctions.Extra+ Numeric.Series Numeric.Sum other-modules: Numeric.SpecFunctions.Internal test-suite tests+ default-language: Haskell2010+ other-extensions: ViewPatterns+ type: exitcode-stdio-1.0 ghc-options: -Wall -threaded if arch(i386)@@ -61,7 +76,7 @@ Tests.Sum build-depends: math-functions,- base >=3 && <5,+ base >=4.5 && <5, deepseq, primitive, vector >= 0.7,
tests/Tests/SpecFunctions.hs view
@@ -14,8 +14,8 @@ import Tests.Helpers import Tests.SpecFunctions.Tables import Numeric.SpecFunctions-import Numeric.MathFunctions.Comparison (within)-+import Numeric.MathFunctions.Comparison (within,relativeError)+import Numeric.MathFunctions.Constants (m_epsilon,m_tiny) tests :: Test tests = testGroup "Special functions"@@ -24,11 +24,11 @@ , testProperty "gamma(1,x) = 1 - exp(-x)" $ incompleteGammaAt1Check , testProperty "0 <= gamma <= 1" $ incompleteGammaInRange , testProperty "0 <= I[B] <= 1" $ incompleteBetaInRange+ , testProperty "invIncompleteGamma = gamma^-1" $ invIGammaIsInverse -- XXX FIXME DISABLED due to failures- -- , testProperty "invIncompleteGamma = gamma^-1" $ invIGammaIsInverse -- , testProperty "invIncompleteBeta = B^-1" $ invIBetaIsInverse- -- , testProperty "gamma - increases" $- -- \s x y -> s > 0 && x > 0 && y > 0 ==> monotonicallyIncreases (incompleteGamma s) x y+ , testProperty "gamma - increases" $+ \(abs -> s) (abs -> x) (abs -> y) -> s > 0 ==> monotonicallyIncreases (incompleteGamma s) x y , testProperty "invErfc = erfc^-1" $ invErfcIsInverse , testProperty "invErf = erf^-1" $ invErfIsInverse -- Unit tests@@ -118,16 +118,22 @@ -- invIncompleteGamma is inverse of incompleteGamma invIGammaIsInverse :: Double -> Double -> Property invIGammaIsInverse (abs -> a) (range01 -> p) =- a > 0 && p > 0 && p < 1 ==> ( counterexample ("a = " ++ show a )- $ counterexample ("p = " ++ show p )- $ counterexample ("x = " ++ show x )- $ counterexample ("p' = " ++ show p')- $ counterexample ("Δp = " ++ show (p - p'))- $ abs (p - p') <= 1e-12- )+ a > m_tiny && p > m_tiny && p < 1 ==>+ ( counterexample ("a = " ++ show a )+ $ counterexample ("p = " ++ show p )+ $ counterexample ("x = " ++ show x )+ $ counterexample ("p' = " ++ show p')+ $ counterexample ("err = " ++ show (relativeError p p'))+ $ counterexample ("pred = " ++ show δ)+ $ relativeError p p' < δ+ ) where x = invIncompleteGamma a p+ f' = exp ( log x * (a-1) - x - logGamma a) p' = incompleteGamma a x+ -- FIXME: 128 is big constant. It should be replaced by something+ -- smaller when #42 is fixed+ δ = (m_epsilon/2) * (256 + 1 * (1 + abs (x * f' / p))) -- invErfc is inverse of erfc invErfcIsInverse :: Double -> Property
− tests/Tests/SpecFunctions_flymake.hs
@@ -1,206 +0,0 @@-{-# LANGUAGE ViewPatterns #-}--- | Tests for Statistics.Math-module Tests.SpecFunctions (- tests- ) where--import qualified Data.Vector as V-import Data.Vector ((!))--import Test.QuickCheck hiding (choose)-import Test.Framework-import Test.Framework.Providers.QuickCheck2--import Tests.Helpers-import Tests.SpecFunctions.Tables-import Numeric.SpecFunctions---tests :: Test-tests = testGroup "Special functions"- [ testProperty "Gamma(x+1) = x*Gamma(x) [logGamma]" $ gammaReccurence logGamma 3e-8- , testProperty "Gamma(x+1) = x*Gamma(x) [logGammaL]" $ gammaReccurence logGammaL 2e-13- , testProperty "gamma(1,x) = 1 - exp(-x)" $ incompleteGammaAt1Check- , testProperty "0 <= gamma <= 1" $ incompleteGammaInRange- , testProperty "0 <= I[B] <= 1" $ incompleteBetaInRange- -- XXX FIXME DISABLED due to failures- -- , testProperty "invIncompleteGamma = gamma^-1" $ invIGammaIsInverse- -- , testProperty "invIncompleteBeta = B^-1" $ invIBetaIsInverse- -- , testProperty "gamma - increases" $- -- \s x y -> s > 0 && x > 0 && y > 0 ==> monotonicallyIncreases (incompleteGamma s) x y- , testProperty "invErfc = erfc^-1" $ invErfcIsInverse- , testProperty "invErf = erf^-1" $ invErfIsInverse- -- Unit tests- , testAssertion "Factorial is expected to be precise at 1e-15 level"- $ and [ eq 1e-15 (factorial (fromIntegral n :: Int))- (fromIntegral (factorial' n))- |n <- [0..170]]- , testAssertion "Log factorial is expected to be precise at 1e-15 level"- $ and [ eq 1e-15 (logFactorial (fromIntegral n :: Int))- (log $ fromIntegral $ factorial' n)- | n <- [2..170]]- , testAssertion "logGamma is expected to be precise at 1e-9 level [integer points]"- $ and [ eq 1e-9 (logGamma (fromIntegral n))- (logFactorial (n-1))- | n <- [3..10000::Int]]- , testAssertion "logGamma is expected to be precise at 1e-9 level [fractional points]"- $ and [ eq 1e-9 (logGamma x) lg | (x,lg) <- tableLogGamma ]- , testAssertion "logGammaL is expected to be precise at 1e-15 level"- $ and [ eq 1e-15 (logGammaL (fromIntegral n))- (logFactorial (n-1))- | n <- [3..10000::Int]]- -- FIXME: Too low!- , testAssertion "logGammaL is expected to be precise at 1e-10 level [fractional points]"- $ and [ eq 1e-10 (logGammaL x) lg | (x,lg) <- tableLogGamma ]- -- FIXME: loss of precision when logBeta p q ≈ 0.- -- Relative error doesn't work properly in this case.- , testAssertion "logBeta is expected to be precise at 1e-6 level"- $ and [ eq 1e-6 (logBeta p q)- (logGammaL p + logGammaL q - logGammaL (p+q))- | p <- [0.1,0.2 .. 0.9] ++ [2 .. 20]- , q <- [0.1,0.2 .. 0.9] ++ [2 .. 20]- ]- , testAssertion "digamma is expected to be precise at 1e-14 [integers]"- $ digammaTestIntegers 1e-14- -- Relative precision is lost when digamma(x) ≈ 0- , testAssertion "digamma is expected to be precise at 1e-12"- $ and [ eq 1e-12 r (digamma x) | (x,r) <- tableDigamma ]- -- FIXME: Why 1e-8? Is it due to poor precision of logBeta?- , testAssertion "incompleteBeta is expected to be precise at 1e-8 level"- $ and [ eq 1e-8 (incompleteBeta p q x) ib | (p,q,x,ib) <- tableIncompleteBeta ]- , testAssertion "incompleteBeta with p > 3000 and q > 3000"- $ and [ eq 1e-11 (incompleteBeta p q x) ib | (x,p,q,ib) <-- [ (0.495, 3001, 3001, 0.2192546757957825068677527085659175689142653854877723)- , (0.501, 3001, 3001, 0.5615652382981522803424365187631195161665429270531389)- , (0.531, 3500, 3200, 0.9209758089734407825580172472327758548870610822321278)- , (0.501, 13500, 13200, 0.0656209987264794057358373443387716674955276089622780)- ]- ]- , testAssertion "choose is expected to precise at 1e-12 level"- $ and [ eq 1e-12 (choose (fromIntegral n) (fromIntegral k)) (fromIntegral $ choose' n k)- | n <- [0..300], k <- [0..n]]- ----------------------------------------------------------------- -- Self tests- , testProperty "Self-test: 0 <= range01 <= 1" $ \x -> let f = range01 x in f <= 1 && f >= 0- ]--------------------------------------------------------------------- QC tests--------------------------------------------------------------------- Γ(x+1) = x·Γ(x)-gammaReccurence :: (Double -> Double) -> Double -> Double -> Property-gammaReccurence logG ε x =- (x > 0 && x < 100) ==> (abs (g2 - g1 - log x) < ε)- where- g1 = logG x- g2 = logG (x+1)---- γ(s,x) is in [0,1] range-incompleteGammaInRange :: Double -> Double -> Property-incompleteGammaInRange (abs -> s) (abs -> x) =- x >= 0 && s > 0 ==> let i = incompleteGamma s x in i >= 0 && i <= 1---- γ(1,x) = 1 - exp(-x)--- Since Γ(1) = 1 normalization doesn't make any difference-incompleteGammaAt1Check :: Double -> Property-incompleteGammaAt1Check (abs -> x) =- x > 0 ==> (incompleteGamma 1 x + exp(-x)) ≈ 1- where- (≈) = eq 1e-13---- invIncompleteGamma is inverse of incompleteGamma-invIGammaIsInverse :: Double -> Double -> Property-invIGammaIsInverse (abs -> a) (range01 -> p) =- a > 0 && p > 0 && p < 1 ==> ( printTestCase ("a = " ++ show a )- $ printTestCase ("p = " ++ show p )- $ printTestCase ("x = " ++ show x )- $ printTestCase ("p' = " ++ show p')- $ printTestCase ("Δp = " ++ show (p - p'))- $ abs (p - p') <= 1e-12- )- where- x = invIncompleteGamma a p- p' = incompleteGamma a x---- invErfc is inverse of erfc-invErfcIsInverse :: Double -> Property-invErfcIsInverse ((*2) . range01 -> p)- = printTestCase ("p = " ++ show p )- $ printTestCase ("x = " ++ show x )- $ printTestCase ("p' = " ++ show p')- $ abs (p - p') <= 1e-14- where- x = invErfc p- p' = erfc x---- invErf is inverse of erf-invErfIsInverse :: Double -> Property-invErfIsInverse a- = printTestCase ("p = " ++ show p )- $ printTestCase ("x = " ++ show x )- $ printTestCase ("p' = " ++ show p')- $ abs (p - p') <= 1e-14- where- x = invErf p- p' = erf x- p | a < 0 = - range01 a- | otherwise = range01 a---- B(s,x) is in [0,1] range-incompleteBetaInRange :: Double -> Double -> Double -> Property-incompleteBetaInRange (abs -> p) (abs -> q) (range01 -> x) =- p > 0 && q > 0 ==> let i = incompleteBeta p q x in i >= 0 && i <= 1---- invIncompleteBeta is inverse of incompleteBeta-invIBetaIsInverse :: Double -> Double -> Double -> Property-invIBetaIsInverse (abs -> p) (abs -> q) (range01 -> x) =- p > 0 && q > 0 ==> ( printTestCase ("p = " ++ show p )- $ printTestCase ("q = " ++ show q )- $ printTestCase ("x = " ++ show x )- $ printTestCase ("x' = " ++ show x')- $ printTestCase ("a = " ++ show a)- $ printTestCase ("err = " ++ (show $ abs $ (x - x') / x))- $ abs (x - x') <= 1e-12- )- where- x' = incompleteBeta p q a- a = invIncompleteBeta p q x---- Table for digamma function:------ Uses equality ψ(n) = H_{n-1} - γ where--- H_{n} = Σ 1/k, k = [1 .. n] - harmonic number--- γ = 0.57721566490153286060 - Euler-Mascheroni number-digammaTestIntegers :: Double -> Bool-digammaTestIntegers eps- = all (uncurry $ eq eps) $ take 3000 digammaInt- where- ok approx exact = approx- -- Harmonic numbers starting from 0- harmN = scanl (\a n -> a + 1/n) 0 [1::Rational .. ]- gam = 0.57721566490153286060- -- Digamma values- digammaInt = zipWith (\i h -> (digamma i, realToFrac h - gam)) [1..] harmN---------------------------------------------------------------------- Unit tests--------------------------------------------------------------------- Lookup table for fact factorial calculation. It has fixed size--- which is bad but it's OK for this particular case-factorial_table :: V.Vector Integer-factorial_table = V.generate 2000 (\n -> product [1..fromIntegral n])---- Exact implementation of factorial-factorial' :: Integer -> Integer-factorial' n = factorial_table ! fromIntegral n---- Exact albeit slow implementation of choose-choose' :: Integer -> Integer -> Integer-choose' n k = factorial' n `div` (factorial' k * factorial' (n-k))---- Truncate double to [0,1]-range01 :: Double -> Double-range01 = abs . (snd :: (Integer, Double) -> Double) . properFraction