diff --git a/Numeric/MathFunctions/Comparison.hs b/Numeric/MathFunctions/Comparison.hs --- a/Numeric/MathFunctions/Comparison.hs +++ b/Numeric/MathFunctions/Comparison.hs @@ -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. diff --git a/Numeric/MathFunctions/Constants.hs b/Numeric/MathFunctions/Constants.hs --- a/Numeric/MathFunctions/Constants.hs +++ b/Numeric/MathFunctions/Constants.hs @@ -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 #-} ---------------------------------------------------------------- diff --git a/Numeric/Polynomial/Chebyshev.hs b/Numeric/Polynomial/Chebyshev.hs --- a/Numeric/Polynomial/Chebyshev.hs +++ b/Numeric/Polynomial/Chebyshev.hs @@ -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 diff --git a/Numeric/RootFinding.hs b/Numeric/RootFinding.hs new file mode 100644 --- /dev/null +++ b/Numeric/RootFinding.hs @@ -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. diff --git a/Numeric/Series.hs b/Numeric/Series.hs new file mode 100644 --- /dev/null +++ b/Numeric/Series.hs @@ -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 diff --git a/Numeric/SpecFunctions.hs b/Numeric/SpecFunctions.hs --- a/Numeric/SpecFunctions.hs +++ b/Numeric/SpecFunctions.hs @@ -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 -- -- --- * Shea, B. (1988) Algorithm AS 239: Chi-squared and incomplete --- gamma integral. /Applied Statistics/ --- 37(3):466–473. +-- * 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. + diff --git a/Numeric/SpecFunctions/Extra.hs b/Numeric/SpecFunctions/Extra.hs --- a/Numeric/SpecFunctions/Extra.hs +++ b/Numeric/SpecFunctions/Extra.hs @@ -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 diff --git a/Numeric/SpecFunctions/Internal.hs b/Numeric/SpecFunctions/Internal.hs --- a/Numeric/SpecFunctions/Internal.hs +++ b/Numeric/SpecFunctions/Internal.hs @@ -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 diff --git a/Numeric/Sum.hs b/Numeric/Sum.hs --- a/Numeric/Sum.hs +++ b/Numeric/Sum.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE BangPatterns, CPP, DeriveDataTypeable, FlexibleContexts, +{-# LANGUAGE BangPatterns, DeriveDataTypeable, FlexibleContexts, MultiParamTypeClasses, TemplateHaskell, TypeFamilies #-} {-# OPTIONS_GHC -fno-warn-name-shadowing #-} -- | diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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 diff --git a/math-functions.cabal b/math-functions.cabal --- a/math-functions.cabal +++ b/math-functions.cabal @@ -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 , @@ -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, diff --git a/tests/Tests/SpecFunctions.hs b/tests/Tests/SpecFunctions.hs --- a/tests/Tests/SpecFunctions.hs +++ b/tests/Tests/SpecFunctions.hs @@ -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 diff --git a/tests/Tests/SpecFunctions_flymake.hs b/tests/Tests/SpecFunctions_flymake.hs deleted file mode 100644 --- a/tests/Tests/SpecFunctions_flymake.hs +++ /dev/null @@ -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