diff --git a/ChangeLog b/ChangeLog deleted file mode 100644 --- a/ChangeLog +++ /dev/null @@ -1,22 +0,0 @@ --*- text -*- - -Changes in 0.1.5 - - * Numeric.Sum: new module adds accurate floating point summation. - -Changes in 0.1.4 - - * logFactorial type is genberalized. It accepts any `Integral' type - - * Evaluation of polynomials using Horner's method where coefficients - are store in lists added - -Changes in 0.1.3 - - * Error function and its inverse added. - - * Digamma function added - - * Evaluation of polynomials using Horner's method added. - - * Crash bug in the inverse incomplete beta fixed. diff --git a/Numeric/SpecFunctions.hs b/Numeric/SpecFunctions.hs --- a/Numeric/SpecFunctions.hs +++ b/Numeric/SpecFunctions.hs @@ -26,6 +26,8 @@ , incompleteBeta , incompleteBeta_ , invIncompleteBeta + -- * Sinc + , sinc -- * Logarithm , log1p , log2 @@ -35,770 +37,12 @@ , stirlingError -- * Combinatorics , choose + , logChoose -- * References -- $references ) where -import Data.Bits ((.&.), (.|.), shiftR) -import Data.Int (Int64) -import qualified Data.Number.Erf as Erf (erfc,erf) -import qualified Data.Vector.Unboxed as U - -import Numeric.Polynomial.Chebyshev (chebyshevBroucke) -import Numeric.Polynomial (evaluateEvenPolynomialL,evaluateOddPolynomialL) -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 - ) -import Text.Printf - - ----------------------------------------------------------------- --- Error function ----------------------------------------------------------------- - --- | Error function. --- --- > erf -∞ = -1 --- > erf 0 = 0 --- > erf +∞ = 1 -erf :: Double -> Double -{-# INLINE erf #-} -erf = Erf.erf - --- | Complementary error function. --- --- > erfc -∞ = 2 --- > erfc 0 = 1 --- > errc +∞ = 0 -erfc :: Double -> Double -{-# INLINE erfc #-} -erfc = Erf.erfc - - --- | Inverse of 'erf'. -invErf :: Double -- ^ /p/ ∈ [-1,1] - -> Double -invErf p = invErfc (1 - p) - --- | Inverse of 'erfc'. -invErfc :: Double -- ^ /p/ ∈ [0,2] - -> Double -invErfc p - | p == 2 = m_neg_inf - | p == 0 = m_pos_inf - | p >0 && p < 2 = if p <= 1 then r else -r - | otherwise = modErr $ "invErfc: p must be in [0,2] got " ++ show p - where - pp = if p <= 1 then p else 2 - p - t = sqrt $ -2 * log( 0.5 * pp) - -- Initial guess - x0 = -0.70711 * ((2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t) - r = loop 0 x0 - -- - loop :: Int -> Double -> Double - loop !j !x - | j >= 2 = x - | otherwise = let err = erfc x - pp - x' = x + err / (1.12837916709551257 * exp(-x * x) - x * err) -- // Halley - in loop (j+1) x' - - - ----------------------------------------------------------------- --- Gamma function ----------------------------------------------------------------- - --- Adapted from http://people.sc.fsu.edu/~burkardt/f_src/asa245/asa245.html - --- | 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). -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) - 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 - --- | 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). -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 - ] - - - --- | Compute the log gamma correction factor 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 -logGammaCorrection :: Double -> Double -logGammaCorrection x - | x < 10 = m_NaN - | x < big = chebyshevBroucke (t * t * 2 - 1) coeffs / x - | otherwise = 1 / (x * 12) - where - big = 94906265.62425156 - t = 10 / x - coeffs = U.fromList [ - 0.1666389480451863247205729650822e+0, - -0.1384948176067563840732986059135e-4, - 0.9810825646924729426157171547487e-8, - -0.1809129475572494194263306266719e-10, - 0.6221098041892605227126015543416e-13, - -0.3399615005417721944303330599666e-15, - 0.2683181998482698748957538846666e-17 - ] - - - --- | Compute the normalized lower incomplete gamma function --- γ(/s/,/x/). Normalization means that --- γ(/s/,∞)=1. Uses Algorithm AS 239 by Shea. -incompleteGamma :: Double -- ^ /s/ ∈ (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 - 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 - -- - 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 - - - --- Adapted from Numerical Recipes §6.2.1 - --- | Inverse incomplete gamma function. It's approximately inverse of --- 'incompleteGamma' for the same /s/. So following equality --- approximately holds: --- --- > invIncompleteGamma s . incompleteGamma s = id -invIncompleteGamma :: Double -- ^ /s/ ∈ (0,∞) - -> Double -- ^ /p/ ∈ [0,1] - -> Double -invIncompleteGamma a p - | a <= 0 = - modErr $ printf "invIncompleteGamma: a must be positive. a=%g p=%g" a p - | p < 0 || p > 1 = - modErr $ printf "invIncompleteGamma: p must be in [0,1] range. a=%g p=%g" a p - | p == 0 = 0 - | p == 1 = 1 / 0 - | otherwise = loop 0 guess - where - -- Solve equation γ(a,x) = p using Halley method - loop :: Int -> Double -> Double - loop i x - | i >= 12 = x' - -- For small s derivative becomes approximately 1/x*exp(-x) and - -- skyrockets for small x. If it happens correct answer is 0. - | isInfinite f' = 0 - | abs dx < eps * x' = x' - | otherwise = loop (i + 1) x' - where - -- Value of γ(a,x) - p - f = incompleteGamma a x - p - -- dγ(a,x)/dx - f' | a > 1 = afac * exp( -(x - a1) + a1 * (log x - lna1)) - | otherwise = exp( -x + a1 * log x - gln) - u = f / f' - -- Halley correction to Newton-Rapson step - corr = u * (a1 / x - 1) - dx = u / (1 - 0.5 * min 1.0 corr) - -- New approximation to x - x' | x < dx = 0.5 * x -- Do not go below 0 - | otherwise = x - dx - -- Calculate inital guess for root - guess - -- - | a > 1 = - let t = sqrt $ -2 * log(if p < 0.5 then p else 1 - p) - x1 = (2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t - x2 = if p < 0.5 then -x1 else x1 - in max 1e-3 (a * (1 - 1/(9*a) - x2 / (3 * sqrt a)) ** 3) - -- For a <= 1 use following approximations: - -- γ(a,1) ≈ 0.253a + 0.12a² - -- - -- γ(a,x) ≈ γ(a,1)·x^a x < 1 - -- γ(a,x) ≈ γ(a,1) + (1 - γ(a,1))(1 - exp(1 - x)) x >= 1 - | otherwise = - let t = 1 - a * (0.253 + a*0.12) - in if p < t - then (p / t) ** (1 / a) - else 1 - log( 1 - (p-t) / (1-t)) - -- Constants - a1 = a - 1 - lna1 = log a1 - afac = exp( a1 * (lna1 - 1) - gln ) - gln = logGamma a - eps = 1e-8 - - - ----------------------------------------------------------------- --- Beta function ----------------------------------------------------------------- - --- | Compute the natural logarithm of the beta function. -logBeta :: Double -> Double -> 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) - | otherwise = logGamma p + logGamma q - logGamma pq - where - p = min a b - q = max a b - ppq = p / pq - 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 - -> Double -- ^ /x/, must lie in [0,1] range - -> Double -incompleteBeta p q = incompleteBeta_ (logBeta p q) p q - --- | 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 -- ^ /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 - | x == 0 || x == 1 = x - | p >= (p+q) * x = incompleteBetaWorker beta p q x - | otherwise = 1 - incompleteBetaWorker beta q p (1 - x) - - --- Approximation of incomplete beta by quandrature. --- --- Note that x =< p/(p+q) -incompleteBetaApprox :: Double -> Double -> Double -> Double -> Double -incompleteBetaApprox beta p q x - | ans > 0 = 1 - ans - | otherwise = -ans - where - -- Constants - p1 = p - 1 - q1 = q - 1 - mu = p / (p + q) - lnmu = log mu - lnmuc = log (1 - mu) - -- Upper limit for integration - xu = max 0 $ min (mu - 10*t) (x - 5*t) - where - t = sqrt $ p*q / ( (p+q) * (p+q) * (p + q + 1) ) - -- Calculate incomplete beta by quadrature - go y w = let t = x + (xu - x) * y - in w * exp( p1 * (log t - lnmu) + q1 * (log(1-t) - lnmuc) ) - s = U.sum $ U.zipWith go coefY coefW - ans = s * (xu - x) * exp( p1 * lnmu + q1 * lnmuc - beta ) - - --- Worker for incomplete beta function. It is separate function to --- avoid confusion with parameter during parameter swapping -incompleteBetaWorker :: Double -> Double -> Double -> Double -> Double -incompleteBetaWorker beta p q x - -- For very large p and q this method becomes very slow so another - -- method is used. - | p > 3000 && q > 3000 = incompleteBetaApprox beta p q x - | otherwise = loop (p+q) (truncate $ q + cx * (p+q)) 1 1 1 - where - -- Constants - eps = 1e-15 - cx = 1 - x - -- Loop - loop !psq (ns :: Int) ai term betain - | done = betain' * exp( p * log x + (q - 1) * log cx - beta) / p - | otherwise = loop psq' (ns - 1) (ai + 1) term' betain' - where - -- New values - term' = term * fact / (p + ai) - betain' = betain + term' - fact | ns > 0 = (q - ai) * x/cx - | ns == 0 = (q - ai) * x - | otherwise = psq * x - -- Iterations are complete - done = db <= eps && db <= eps*betain' where db = abs term' - psq' = if ns < 0 then psq + 1 else psq - - - --- | 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] - -> Double -invIncompleteBeta p q a - | p <= 0 || q <= 0 = - modErr $ printf "invIncompleteBeta p <= 0 || q <= 0. p=%g q=%g a=%g" p q a - | a < 0 || a > 1 = - modErr $ printf "invIncompleteBeta x must be in [0,1]. p=%g q=%g a=%g" p q a - | a == 0 || a == 1 = a - | a > 0.5 = 1 - invIncompleteBetaWorker (logBeta p q) q p (1 - a) - | otherwise = invIncompleteBetaWorker (logBeta p q) p q a - - -invIncompleteBetaWorker :: Double -> Double -> Double -> Double -> Double --- NOTE: p <= 0.5. -invIncompleteBetaWorker beta a b p = loop (0::Int) guess - where - a1 = a - 1 - b1 = b - 1 - -- Solve equation using Halley method - loop !i !x - -- We cannot continue at this point so we simply return `x' - | x == 0 || x == 1 = x - -- When derivative becomes infinite we cannot continue - -- iterations. It can only happen in vicinity of 0 or 1. It's - -- hardly possible to get good answer in such circumstances but - -- `x' is already reasonable. - | isInfinite f' = x - -- Iterations limit reached. Most of the time solution will - -- converge to answer because of discreteness of Double. But - -- solution have good precision already. - | i >= 10 = x - -- Solution converges - | abs dx <= 16 * m_epsilon * x = x' - | otherwise = loop (i+1) x' - where - -- Calculate Halley step. - f = incompleteBeta_ beta a b x - p - f' = exp $ a1 * log x + b1 * log (1 - x) - beta - u = f / f' - dx = u / (1 - 0.5 * min 1 (u * (a1 / x - b1 / (1 - x)))) - -- 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 - 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 - - - - ----------------------------------------------------------------- --- Logarithm ----------------------------------------------------------------- - --- | Compute the natural logarithm of 1 + @x@. This is accurate even --- for values of @x@ near zero, where use of @log(1+x)@ would lose --- precision. -log1p :: Double -> Double -log1p x - | x == 0 = 0 - | x == -1 = m_neg_inf - | x < -1 = m_NaN - | x' < m_epsilon * 0.5 = x - | (x >= 0 && x < 1e-8) || (x >= -1e-9 && x < 0) - = x * (1 - x * 0.5) - | x' < 0.375 = x * (1 - x * chebyshevBroucke (x / 0.375) coeffs) - | otherwise = log (1 + x) - where - x' = abs x - coeffs = U.fromList [ - 0.10378693562743769800686267719098e+1, - -0.13364301504908918098766041553133e+0, - 0.19408249135520563357926199374750e-1, - -0.30107551127535777690376537776592e-2, - 0.48694614797154850090456366509137e-3, - -0.81054881893175356066809943008622e-4, - 0.13778847799559524782938251496059e-4, - -0.23802210894358970251369992914935e-5, - 0.41640416213865183476391859901989e-6, - -0.73595828378075994984266837031998e-7, - 0.13117611876241674949152294345011e-7, - -0.23546709317742425136696092330175e-8, - 0.42522773276034997775638052962567e-9, - -0.77190894134840796826108107493300e-10, - 0.14075746481359069909215356472191e-10, - -0.25769072058024680627537078627584e-11, - 0.47342406666294421849154395005938e-12, - -0.87249012674742641745301263292675e-13, - 0.16124614902740551465739833119115e-13, - -0.29875652015665773006710792416815e-14, - 0.55480701209082887983041321697279e-15, - -0.10324619158271569595141333961932e-15 - ] - - --- | /O(log n)/ Compute the logarithm in base 2 of the given value. -log2 :: Int -> Int -log2 v0 - | v0 <= 0 = modErr $ "log2: negative input, got " ++ show v0 - | otherwise = go 5 0 v0 - where - go !i !r !v | i == -1 = r - | v .&. b i /= 0 = let si = U.unsafeIndex sv i - 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] - !sv = U.fromList [1,2,4,8,16,32] - - ----------------------------------------------------------------- --- Factorial ----------------------------------------------------------------- - --- | Compute the factorial function /n/!. Returns +∞ if the --- input is above 170 (above which the result cannot be represented by --- a 64-bit 'Double'). -factorial :: Int -> Double -factorial n - | n < 0 = error "Numeric.SpecFunctions.factorial: negative input" - | n <= 1 = 1 - | n <= 170 = U.product $ U.map fromIntegral $ U.enumFromTo 2 n - | otherwise = m_pos_inf - --- | Compute the natural logarithm of the factorial function. Gives --- 16 decimal digits of precision. -logFactorial :: Integral a => a -> Double -logFactorial n - | n < 0 = error "Numeric.SpecFunctions.logFactorial: negative input" - | n <= 14 = log $ factorial $ fromIntegral n - | otherwise = (x - 0.5) * log x - x + 9.1893853320467e-1 + z / x - where x = fromIntegral n + 1 - y = 1 / (x * x) - z = ((-(5.95238095238e-4 * y) + 7.936500793651e-4) * y - - 2.7777777777778e-3) * y + 8.3333333333333e-2 -{-# SPECIALIZE logFactorial :: Int -> Double #-} - --- | Calculate the error term of the Stirling approximation. This is --- only defined for non-negative values. --- --- > stirlingError @n@ = @log(n!) - log(sqrt(2*pi*n)*(n/e)^n) -stirlingError :: Double -> Double -stirlingError n - | n <= 15.0 = case properFraction (n+n) of - (i,0) -> sfe `U.unsafeIndex` i - _ -> logGamma (n+1.0) - (n+0.5) * log n + n - - m_ln_sqrt_2_pi - | n > 500 = evaluateOddPolynomialL (1/n) [s0,-s1] - | n > 80 = evaluateOddPolynomialL (1/n) [s0,-s1,s2] - | n > 35 = evaluateOddPolynomialL (1/n) [s0,-s1,s2,-s3] - | otherwise = evaluateOddPolynomialL (1/n) [s0,-s1,s2,-s3,s4] - where - s0 = 0.083333333333333333333 -- 1/12 - s1 = 0.00277777777777777777778 -- 1/360 - s2 = 0.00079365079365079365079365 -- 1/1260 - s3 = 0.000595238095238095238095238 -- 1/1680 - s4 = 0.0008417508417508417508417508 -- 1/1188 - sfe = U.fromList [ 0.0, - 0.1534264097200273452913848, 0.0810614667953272582196702, - 0.0548141210519176538961390, 0.0413406959554092940938221, - 0.03316287351993628748511048, 0.02767792568499833914878929, - 0.02374616365629749597132920, 0.02079067210376509311152277, - 0.01848845053267318523077934, 0.01664469118982119216319487, - 0.01513497322191737887351255, 0.01387612882307074799874573, - 0.01281046524292022692424986, 0.01189670994589177009505572, - 0.01110455975820691732662991, 0.010411265261972096497478567, - 0.009799416126158803298389475, 0.009255462182712732917728637, - 0.008768700134139385462952823, 0.008330563433362871256469318, - 0.007934114564314020547248100, 0.007573675487951840794972024, - 0.007244554301320383179543912, 0.006942840107209529865664152, - 0.006665247032707682442354394, 0.006408994188004207068439631, - 0.006171712263039457647532867, 0.005951370112758847735624416, - 0.005746216513010115682023589, 0.005554733551962801371038690 ] - - ----------------------------------------------------------------- --- Combinatorics ----------------------------------------------------------------- - --- | Quickly compute the natural logarithm of /n/ @`choose`@ /k/, with --- no checking. -logChooseFast :: Double -> Double -> Double -logChooseFast n k = -log (n + 1) - logBeta (n - k + 1) (k + 1) - --- | Compute the binomial coefficient /n/ @\``choose`\`@ /k/. For --- values of /k/ > 30, this uses an approximation for performance --- reasons. The approximation is accurate to 12 decimal places in the --- worst case --- --- Example: --- --- > 7 `choose` 3 == 35 -choose :: Int -> Int -> Double -n `choose` k - | k > n = 0 - | k' < 50 = U.foldl' go 1 . U.enumFromTo 1 $ k' - | approx < max64 = fromIntegral . round64 $ approx - | otherwise = approx - where - k' = min k (n-k) - approx = exp $ logChooseFast (fromIntegral n) (fromIntegral k') - -- Less numerically stable: - -- exp $ lg (n+1) - lg (k+1) - lg (n-k+1) - -- where lg = logGamma . fromIntegral - go a i = a * (nk + j) / j - where j = fromIntegral i :: Double - nk = fromIntegral (n - k') - 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. -digamma :: Double -> Double -digamma x - | isNaN x || isInfinite x = m_NaN - -- FIXME: - -- This is ugly. We are testing here that number is in fact - -- integer. It's somewhat tricky question to answer. When ε for - -- given number becomes 1 or greater every number is represents - -- an integer. We also must make sure that excess precision - -- won't bite us. - | x <= 0 && fromIntegral (truncate x :: Int64) == x = m_neg_inf - -- Jeffery's reflection formula - | x < 0 = digamma (1 - x) + pi / tan (negate pi * x) - | x <= 1e-6 = - γ - 1/x + trigamma1 * x - | x' < c = r - -- De Moivre's expansion - | otherwise = let s = 1/x' - in evaluateEvenPolynomialL s - [ r + log x' - 0.5 * s - , - 1/12 - , 1/120 - , - 1/252 - , 1/240 - , - 1/132 - , 391/32760 - ] - where - γ = m_eulerMascheroni - c = 12 - -- Reduce to digamma (x + n) where (x + n) >= c - (r, x') = reduce 0 x - where - reduce !s y - | y < c = reduce (s - 1 / y) (y + 1) - | otherwise = (s, y) - - - ----------------------------------------------------------------- --- Constants ----------------------------------------------------------------- - --- Coefficients for 18-point Gauss-Legendre integration. They are --- used in implementation of incomplete gamma and beta functions. -coefW,coefY :: U.Vector Double -coefW = U.fromList [ 0.0055657196642445571, 0.012915947284065419, 0.020181515297735382 - , 0.027298621498568734, 0.034213810770299537, 0.040875750923643261 - , 0.047235083490265582, 0.053244713977759692, 0.058860144245324798 - , 0.064039797355015485, 0.068745323835736408, 0.072941885005653087 - , 0.076598410645870640, 0.079687828912071670, 0.082187266704339706 - , 0.084078218979661945, 0.085346685739338721, 0.085983275670394821 - ] -coefY = U.fromList [ 0.0021695375159141994, 0.011413521097787704, 0.027972308950302116 - , 0.051727015600492421, 0.082502225484340941, 0.12007019910960293 - , 0.16415283300752470, 0.21442376986779355, 0.27051082840644336 - , 0.33199876341447887, 0.39843234186401943, 0.46931971407375483 - , 0.54413605556657973, 0.62232745288031077, 0.70331500465597174 - , 0.78649910768313447, 0.87126389619061517, 0.95698180152629142 - ] -{-# NOINLINE coefW #-} -{-# NOINLINE coefY #-} - -trigamma1 :: Double -trigamma1 = 1.6449340668482264365 -- pi**2 / 6 - -modErr :: String -> a -modErr msg = error $ "Numeric.SpecFunctions." ++ msg - - +import Numeric.SpecFunctions.Internal -- $references -- diff --git a/Numeric/SpecFunctions/Extra.hs b/Numeric/SpecFunctions/Extra.hs --- a/Numeric/SpecFunctions/Extra.hs +++ b/Numeric/SpecFunctions/Extra.hs @@ -9,10 +9,13 @@ -- -- Less common mathematical functions. module Numeric.SpecFunctions.Extra ( - bd0 + bd0 + , chooseExact + , logChooseFast ) where import Numeric.MathFunctions.Constants (m_NaN) +import Numeric.SpecFunctions.Internal (chooseExact,logChooseFast) -- | Evaluate the deviance term @x log(x/np) + np - x@. bd0 :: Double -- ^ @x@ diff --git a/Numeric/SpecFunctions/Internal.hs b/Numeric/SpecFunctions/Internal.hs new file mode 100644 --- /dev/null +++ b/Numeric/SpecFunctions/Internal.hs @@ -0,0 +1,810 @@ +{-# LANGUAGE BangPatterns, ScopedTypeVariables #-} +-- | +-- Module : Numeric.SpecFunctions.Internal +-- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Internal module with implementation of special functions. +module Numeric.SpecFunctions.Internal where + +import Data.Bits ((.&.), (.|.), shiftR) +import Data.Int (Int64) +import qualified Data.Number.Erf as Erf (erfc,erf) +import qualified Data.Vector.Unboxed as U + +import Numeric.Polynomial.Chebyshev (chebyshevBroucke) +import Numeric.Polynomial (evaluateEvenPolynomialL,evaluateOddPolynomialL) +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 + ) +import Text.Printf + + +---------------------------------------------------------------- +-- Error function +---------------------------------------------------------------- + +-- | Error function. +-- +-- > erf -∞ = -1 +-- > erf 0 = 0 +-- > erf +∞ = 1 +erf :: Double -> Double +{-# INLINE erf #-} +erf = Erf.erf + +-- | Complementary error function. +-- +-- > erfc -∞ = 2 +-- > erfc 0 = 1 +-- > errc +∞ = 0 +erfc :: Double -> Double +{-# INLINE erfc #-} +erfc = Erf.erfc + + +-- | Inverse of 'erf'. +invErf :: Double -- ^ /p/ ∈ [-1,1] + -> Double +invErf p = invErfc (1 - p) + +-- | Inverse of 'erfc'. +invErfc :: Double -- ^ /p/ ∈ [0,2] + -> Double +invErfc p + | p == 2 = m_neg_inf + | p == 0 = m_pos_inf + | p >0 && p < 2 = if p <= 1 then r else -r + | otherwise = modErr $ "invErfc: p must be in [0,2] got " ++ show p + where + pp = if p <= 1 then p else 2 - p + t = sqrt $ -2 * log( 0.5 * pp) + -- Initial guess + x0 = -0.70711 * ((2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t) + r = loop 0 x0 + -- + loop :: Int -> Double -> Double + loop !j !x + | j >= 2 = x + | otherwise = let err = erfc x - pp + x' = x + err / (1.12837916709551257 * exp(-x * x) - x * err) -- // Halley + in loop (j+1) x' + + + +---------------------------------------------------------------- +-- Gamma function +---------------------------------------------------------------- + +-- Adapted from http://people.sc.fsu.edu/~burkardt/f_src/asa245/asa245.html + +-- | 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). +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) + 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 + +-- | 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). +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 + ] + + + +-- | Compute the log gamma correction factor 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 +logGammaCorrection :: Double -> Double +logGammaCorrection x + | x < 10 = m_NaN + | x < big = chebyshevBroucke (t * t * 2 - 1) coeffs / x + | otherwise = 1 / (x * 12) + where + big = 94906265.62425156 + t = 10 / x + coeffs = U.fromList [ + 0.1666389480451863247205729650822e+0, + -0.1384948176067563840732986059135e-4, + 0.9810825646924729426157171547487e-8, + -0.1809129475572494194263306266719e-10, + 0.6221098041892605227126015543416e-13, + -0.3399615005417721944303330599666e-15, + 0.2683181998482698748957538846666e-17 + ] + + + +-- | Compute the normalized lower incomplete gamma function +-- γ(/s/,/x/). Normalization means that +-- γ(/s/,∞)=1. Uses Algorithm AS 239 by Shea. +incompleteGamma :: Double -- ^ /s/ ∈ (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 + 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 + -- + 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 + + + +-- Adapted from Numerical Recipes §6.2.1 + +-- | Inverse incomplete gamma function. It's approximately inverse of +-- 'incompleteGamma' for the same /s/. So following equality +-- approximately holds: +-- +-- > invIncompleteGamma s . incompleteGamma s = id +invIncompleteGamma :: Double -- ^ /s/ ∈ (0,∞) + -> Double -- ^ /p/ ∈ [0,1] + -> Double +invIncompleteGamma a p + | a <= 0 = + modErr $ printf "invIncompleteGamma: a must be positive. a=%g p=%g" a p + | p < 0 || p > 1 = + modErr $ printf "invIncompleteGamma: p must be in [0,1] range. a=%g p=%g" a p + | p == 0 = 0 + | p == 1 = 1 / 0 + | otherwise = loop 0 guess + where + -- Solve equation γ(a,x) = p using Halley method + loop :: Int -> Double -> Double + loop i x + | i >= 12 = x' + -- For small s derivative becomes approximately 1/x*exp(-x) and + -- skyrockets for small x. If it happens correct answer is 0. + | isInfinite f' = 0 + | abs dx < eps * x' = x' + | otherwise = loop (i + 1) x' + where + -- Value of γ(a,x) - p + f = incompleteGamma a x - p + -- dγ(a,x)/dx + f' | a > 1 = afac * exp( -(x - a1) + a1 * (log x - lna1)) + | otherwise = exp( -x + a1 * log x - gln) + u = f / f' + -- Halley correction to Newton-Rapson step + corr = u * (a1 / x - 1) + dx = u / (1 - 0.5 * min 1.0 corr) + -- New approximation to x + x' | x < dx = 0.5 * x -- Do not go below 0 + | otherwise = x - dx + -- Calculate inital guess for root + guess + -- + | a > 1 = + let t = sqrt $ -2 * log(if p < 0.5 then p else 1 - p) + x1 = (2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t + x2 = if p < 0.5 then -x1 else x1 + in max 1e-3 (a * (1 - 1/(9*a) - x2 / (3 * sqrt a)) ** 3) + -- For a <= 1 use following approximations: + -- γ(a,1) ≈ 0.253a + 0.12a² + -- + -- γ(a,x) ≈ γ(a,1)·x^a x < 1 + -- γ(a,x) ≈ γ(a,1) + (1 - γ(a,1))(1 - exp(1 - x)) x >= 1 + | otherwise = + let t = 1 - a * (0.253 + a*0.12) + in if p < t + then (p / t) ** (1 / a) + else 1 - log( 1 - (p-t) / (1-t)) + -- Constants + a1 = a - 1 + lna1 = log a1 + afac = exp( a1 * (lna1 - 1) - gln ) + gln = logGamma a + eps = 1e-8 + + + +---------------------------------------------------------------- +-- Beta function +---------------------------------------------------------------- + +-- | Compute the natural logarithm of the beta function. +logBeta :: Double -> Double -> 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) + | otherwise = logGamma p + logGamma q - logGamma pq + where + p = min a b + q = max a b + ppq = p / pq + 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 + -> Double -- ^ /x/, must lie in [0,1] range + -> Double +incompleteBeta p q = incompleteBeta_ (logBeta p q) p q + +-- | 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 -- ^ /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 + | x == 0 || x == 1 = x + | p >= (p+q) * x = incompleteBetaWorker beta p q x + | otherwise = 1 - incompleteBetaWorker beta q p (1 - x) + + +-- Approximation of incomplete beta by quandrature. +-- +-- Note that x =< p/(p+q) +incompleteBetaApprox :: Double -> Double -> Double -> Double -> Double +incompleteBetaApprox beta p q x + | ans > 0 = 1 - ans + | otherwise = -ans + where + -- Constants + p1 = p - 1 + q1 = q - 1 + mu = p / (p + q) + lnmu = log mu + lnmuc = log (1 - mu) + -- Upper limit for integration + xu = max 0 $ min (mu - 10*t) (x - 5*t) + where + t = sqrt $ p*q / ( (p+q) * (p+q) * (p + q + 1) ) + -- Calculate incomplete beta by quadrature + go y w = let t = x + (xu - x) * y + in w * exp( p1 * (log t - lnmu) + q1 * (log(1-t) - lnmuc) ) + s = U.sum $ U.zipWith go coefY coefW + ans = s * (xu - x) * exp( p1 * lnmu + q1 * lnmuc - beta ) + + +-- Worker for incomplete beta function. It is separate function to +-- avoid confusion with parameter during parameter swapping +incompleteBetaWorker :: Double -> Double -> Double -> Double -> Double +incompleteBetaWorker beta p q x + -- For very large p and q this method becomes very slow so another + -- method is used. + | p > 3000 && q > 3000 = incompleteBetaApprox beta p q x + | otherwise = loop (p+q) (truncate $ q + cx * (p+q)) 1 1 1 + where + -- Constants + eps = 1e-15 + cx = 1 - x + -- Loop + loop !psq (ns :: Int) ai term betain + | done = betain' * exp( p * log x + (q - 1) * log cx - beta) / p + | otherwise = loop psq' (ns - 1) (ai + 1) term' betain' + where + -- New values + term' = term * fact / (p + ai) + betain' = betain + term' + fact | ns > 0 = (q - ai) * x/cx + | ns == 0 = (q - ai) * x + | otherwise = psq * x + -- Iterations are complete + done = db <= eps && db <= eps*betain' where db = abs term' + psq' = if ns < 0 then psq + 1 else psq + + + +-- | 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] + -> Double +invIncompleteBeta p q a + | p <= 0 || q <= 0 = + modErr $ printf "invIncompleteBeta p <= 0 || q <= 0. p=%g q=%g a=%g" p q a + | a < 0 || a > 1 = + modErr $ printf "invIncompleteBeta x must be in [0,1]. p=%g q=%g a=%g" p q a + | a == 0 || a == 1 = a + | a > 0.5 = 1 - invIncompleteBetaWorker (logBeta p q) q p (1 - a) + | otherwise = invIncompleteBetaWorker (logBeta p q) p q a + + +invIncompleteBetaWorker :: Double -> Double -> Double -> Double -> Double +-- NOTE: p <= 0.5. +invIncompleteBetaWorker beta a b p = loop (0::Int) guess + where + a1 = a - 1 + b1 = b - 1 + -- Solve equation using Halley method + loop !i !x + -- We cannot continue at this point so we simply return `x' + | x == 0 || x == 1 = x + -- When derivative becomes infinite we cannot continue + -- iterations. It can only happen in vicinity of 0 or 1. It's + -- hardly possible to get good answer in such circumstances but + -- `x' is already reasonable. + | isInfinite f' = x + -- Iterations limit reached. Most of the time solution will + -- converge to answer because of discreteness of Double. But + -- solution have good precision already. + | i >= 10 = x + -- Solution converges + | abs dx <= 16 * m_epsilon * x = x' + | otherwise = loop (i+1) x' + where + -- Calculate Halley step. + f = incompleteBeta_ beta a b x - p + f' = exp $ a1 * log x + b1 * log (1 - x) - beta + u = f / f' + dx = u / (1 - 0.5 * min 1 (u * (a1 / x - b1 / (1 - x)))) + -- 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 + 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 + + + +---------------------------------------------------------------- +-- Sinc function +---------------------------------------------------------------- + +-- | Compute sinc function @sin(x)\/x@ +sinc :: Double -> Double +sinc x + | ax < eps_0 = 1 + | ax < eps_2 = 1 - x2/6 + | ax < eps_4 = 1 - x2/6 + x2*x2/120 + | otherwise = sin x / x + where + ax = abs x + x2 = x*x + -- For explanation of choice see `doc/sinc.hs' + eps_0 = 1.8250120749944284e-8 -- sqrt (6ε/4) + eps_2 = 1.4284346431400855e-4 -- (30ε)**(1/4) / 2 + eps_4 = 4.043633626430947e-3 -- (1206ε)**(1/6) / 2 + + +---------------------------------------------------------------- +-- Logarithm +---------------------------------------------------------------- + +-- | Compute the natural logarithm of 1 + @x@. This is accurate even +-- for values of @x@ near zero, where use of @log(1+x)@ would lose +-- precision. +log1p :: Double -> Double +log1p x + | x == 0 = 0 + | x == -1 = m_neg_inf + | x < -1 = m_NaN + | x' < m_epsilon * 0.5 = x + | (x >= 0 && x < 1e-8) || (x >= -1e-9 && x < 0) + = x * (1 - x * 0.5) + | x' < 0.375 = x * (1 - x * chebyshevBroucke (x / 0.375) coeffs) + | otherwise = log (1 + x) + where + x' = abs x + coeffs = U.fromList [ + 0.10378693562743769800686267719098e+1, + -0.13364301504908918098766041553133e+0, + 0.19408249135520563357926199374750e-1, + -0.30107551127535777690376537776592e-2, + 0.48694614797154850090456366509137e-3, + -0.81054881893175356066809943008622e-4, + 0.13778847799559524782938251496059e-4, + -0.23802210894358970251369992914935e-5, + 0.41640416213865183476391859901989e-6, + -0.73595828378075994984266837031998e-7, + 0.13117611876241674949152294345011e-7, + -0.23546709317742425136696092330175e-8, + 0.42522773276034997775638052962567e-9, + -0.77190894134840796826108107493300e-10, + 0.14075746481359069909215356472191e-10, + -0.25769072058024680627537078627584e-11, + 0.47342406666294421849154395005938e-12, + -0.87249012674742641745301263292675e-13, + 0.16124614902740551465739833119115e-13, + -0.29875652015665773006710792416815e-14, + 0.55480701209082887983041321697279e-15, + -0.10324619158271569595141333961932e-15 + ] + + +-- | /O(log n)/ Compute the logarithm in base 2 of the given value. +log2 :: Int -> Int +log2 v0 + | v0 <= 0 = modErr $ "log2: nonpositive input, got " ++ show v0 + | otherwise = go 5 0 v0 + where + go !i !r !v | i == -1 = r + | v .&. b i /= 0 = let si = U.unsafeIndex sv i + 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] + !sv = U.fromList [1,2,4,8,16,32] + + +---------------------------------------------------------------- +-- Factorial +---------------------------------------------------------------- + +-- | Compute the factorial function /n/!. Returns +∞ if the +-- input is above 170 (above which the result cannot be represented by +-- a 64-bit 'Double'). +factorial :: Int -> Double +factorial n + | n < 0 = error "Numeric.SpecFunctions.factorial: negative input" + | n <= 1 = 1 + | n <= 170 = U.product $ U.map fromIntegral $ U.enumFromTo 2 n + | otherwise = m_pos_inf + +-- | Compute the natural logarithm of the factorial function. Gives +-- 16 decimal digits of precision. +logFactorial :: Integral a => a -> Double +logFactorial n + | n < 0 = error "Numeric.SpecFunctions.logFactorial: negative input" + | n <= 14 = log $ factorial $ fromIntegral n + | otherwise = (x - 0.5) * log x - x + 9.1893853320467e-1 + z / x + where x = fromIntegral n + 1 + y = 1 / (x * x) + z = ((-(5.95238095238e-4 * y) + 7.936500793651e-4) * y - + 2.7777777777778e-3) * y + 8.3333333333333e-2 +{-# SPECIALIZE logFactorial :: Int -> Double #-} + +-- | Calculate the error term of the Stirling approximation. This is +-- only defined for non-negative values. +-- +-- > stirlingError @n@ = @log(n!) - log(sqrt(2*pi*n)*(n/e)^n) +stirlingError :: Double -> Double +stirlingError n + | n <= 15.0 = case properFraction (n+n) of + (i,0) -> sfe `U.unsafeIndex` i + _ -> logGamma (n+1.0) - (n+0.5) * log n + n - + m_ln_sqrt_2_pi + | n > 500 = evaluateOddPolynomialL (1/n) [s0,-s1] + | n > 80 = evaluateOddPolynomialL (1/n) [s0,-s1,s2] + | n > 35 = evaluateOddPolynomialL (1/n) [s0,-s1,s2,-s3] + | otherwise = evaluateOddPolynomialL (1/n) [s0,-s1,s2,-s3,s4] + where + s0 = 0.083333333333333333333 -- 1/12 + s1 = 0.00277777777777777777778 -- 1/360 + s2 = 0.00079365079365079365079365 -- 1/1260 + s3 = 0.000595238095238095238095238 -- 1/1680 + s4 = 0.0008417508417508417508417508 -- 1/1188 + sfe = U.fromList [ 0.0, + 0.1534264097200273452913848, 0.0810614667953272582196702, + 0.0548141210519176538961390, 0.0413406959554092940938221, + 0.03316287351993628748511048, 0.02767792568499833914878929, + 0.02374616365629749597132920, 0.02079067210376509311152277, + 0.01848845053267318523077934, 0.01664469118982119216319487, + 0.01513497322191737887351255, 0.01387612882307074799874573, + 0.01281046524292022692424986, 0.01189670994589177009505572, + 0.01110455975820691732662991, 0.010411265261972096497478567, + 0.009799416126158803298389475, 0.009255462182712732917728637, + 0.008768700134139385462952823, 0.008330563433362871256469318, + 0.007934114564314020547248100, 0.007573675487951840794972024, + 0.007244554301320383179543912, 0.006942840107209529865664152, + 0.006665247032707682442354394, 0.006408994188004207068439631, + 0.006171712263039457647532867, 0.005951370112758847735624416, + 0.005746216513010115682023589, 0.005554733551962801371038690 ] + + +---------------------------------------------------------------- +-- Combinatorics +---------------------------------------------------------------- + +-- | +-- Quickly compute the natural logarithm of /n/ @`choose`@ /k/, with +-- no checking. +-- +-- Less numerically stable: +-- +-- > exp $ lg (n+1) - lg (k+1) - lg (n-k+1) +-- > where lg = logGamma . fromIntegral +logChooseFast :: Double -> Double -> Double +logChooseFast n k = -log (n + 1) - logBeta (n - k + 1) (k + 1) + +-- | Calculate binomial coefficient using exact formula +chooseExact :: Int -> Int -> Double +n `chooseExact` k + = U.foldl' go 1 $ U.enumFromTo 1 k + where + go a i = a * (nk + j) / j + where j = fromIntegral i :: Double + nk = fromIntegral (n - k) + +-- | Compute logarithm of the binomial coefficient. +logChoose :: Int -> Int -> Double +n `logChoose` k + | k > n = (-1) / 0 + -- For very large N exact algorithm overflows double so we + -- switch to beta-function based one + | k' < 50 && (n < 20000000) = log $ chooseExact n k' + | otherwise = logChooseFast (fromIntegral n) (fromIntegral k) + where + k' = min k (n-k) + +-- | Compute the binomial coefficient /n/ @\``choose`\`@ /k/. For +-- values of /k/ > 50, this uses an approximation for performance +-- reasons. The approximation is accurate to 12 decimal places in the +-- worst case +-- +-- Example: +-- +-- > 7 `choose` 3 == 35 +choose :: Int -> Int -> Double +n `choose` k + | k > n = 0 + | k' < 50 = chooseExact n k' + | approx < max64 = fromIntegral . round64 $ approx + | otherwise = approx + where + k' = min k (n-k) + approx = exp $ logChooseFast (fromIntegral n) (fromIntegral k') + 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. +digamma :: Double -> Double +digamma x + | isNaN x || isInfinite x = m_NaN + -- FIXME: + -- This is ugly. We are testing here that number is in fact + -- integer. It's somewhat tricky question to answer. When ε for + -- given number becomes 1 or greater every number is represents + -- an integer. We also must make sure that excess precision + -- won't bite us. + | x <= 0 && fromIntegral (truncate x :: Int64) == x = m_neg_inf + -- Jeffery's reflection formula + | x < 0 = digamma (1 - x) + pi / tan (negate pi * x) + | x <= 1e-6 = - γ - 1/x + trigamma1 * x + | x' < c = r + -- De Moivre's expansion + | otherwise = let s = 1/x' + in evaluateEvenPolynomialL s + [ r + log x' - 0.5 * s + , - 1/12 + , 1/120 + , - 1/252 + , 1/240 + , - 1/132 + , 391/32760 + ] + where + γ = m_eulerMascheroni + c = 12 + -- Reduce to digamma (x + n) where (x + n) >= c + (r, x') = reduce 0 x + where + reduce !s y + | y < c = reduce (s - 1 / y) (y + 1) + | otherwise = (s, y) + + + +---------------------------------------------------------------- +-- Constants +---------------------------------------------------------------- + +-- Coefficients for 18-point Gauss-Legendre integration. They are +-- used in implementation of incomplete gamma and beta functions. +coefW,coefY :: U.Vector Double +coefW = U.fromList [ 0.0055657196642445571, 0.012915947284065419, 0.020181515297735382 + , 0.027298621498568734, 0.034213810770299537, 0.040875750923643261 + , 0.047235083490265582, 0.053244713977759692, 0.058860144245324798 + , 0.064039797355015485, 0.068745323835736408, 0.072941885005653087 + , 0.076598410645870640, 0.079687828912071670, 0.082187266704339706 + , 0.084078218979661945, 0.085346685739338721, 0.085983275670394821 + ] +coefY = U.fromList [ 0.0021695375159141994, 0.011413521097787704, 0.027972308950302116 + , 0.051727015600492421, 0.082502225484340941, 0.12007019910960293 + , 0.16415283300752470, 0.21442376986779355, 0.27051082840644336 + , 0.33199876341447887, 0.39843234186401943, 0.46931971407375483 + , 0.54413605556657973, 0.62232745288031077, 0.70331500465597174 + , 0.78649910768313447, 0.87126389619061517, 0.95698180152629142 + ] +{-# NOINLINE coefW #-} +{-# NOINLINE coefY #-} + +trigamma1 :: Double +trigamma1 = 1.6449340668482264365 -- pi**2 / 6 + +modErr :: String -> a +modErr msg = error $ "Numeric.SpecFunctions." ++ msg diff --git a/benchmark/bench.hs b/benchmark/bench.hs --- a/benchmark/bench.hs +++ b/benchmark/bench.hs @@ -76,6 +76,11 @@ , 100 ] ] + , bgroup "sinc" $ + bench "sin" (nf sin (0.55 :: Double)) + : [ bench (show x) $ nf sinc x + | x <- [0, 1e-6, 1e-3, 0.5] + ] , bgroup "poly" $ [ bench ("vector_"++show (U.length coefs)) $ nf (\x -> evaluatePolynomial x coefs) (1 :: Double) | coefs <- coef_list ] diff --git a/changelog.md b/changelog.md new file mode 100644 --- /dev/null +++ b/changelog.md @@ -0,0 +1,32 @@ +Changes in 0.1.6.0 + + * `logChoose` added for calculation of logarithm of binomial coefficient + + * `chooseExact` and `logChooseFast` added + + * `sinc` added + +Changes in 0.1.5.3 + + * Fix for test suite on 32bit platform + +Changes in 0.1.5 + + * Numeric.Sum: new module adds accurate floating point summation. + +Changes in 0.1.4 + + * logFactorial type is genberalized. It accepts any `Integral` type + + * Evaluation of polynomials using Horner's method where coefficients + are store in lists added + +Changes in 0.1.3 + + * Error function and its inverse added. + + * Digamma function added + + * Evaluation of polynomials using Horner's method added. + + * Crash bug in the inverse incomplete beta fixed. diff --git a/doc/sinc.hs b/doc/sinc.hs new file mode 100644 --- /dev/null +++ b/doc/sinc.hs @@ -0,0 +1,33 @@ +-- Description of choice of approximation boundaries in sinc function +module Sinc where + +import Numeric.MathFunctions.Constants (m_epsilon) + + +-- Approximations for sinc up to 6th order and "exact" implementation +f2,f4,f6,f :: Double -> Double +f2 x = 1 - x*x/6 +f4 x = 1 - x*x/6 + x*x*x*x/120 +f6 x = 1 - x*x/6 + x*x*x*x/120 - x*x*x*x*x*x/5040 +f x = sin x / x + +-- When next term becomes so small that (1-e)==1 we can neglect it: +e0,e2,e4 :: Double +e0 = sqrt (6 * m_epsilon / 4) +e2 = (30 * m_epsilon) ** (1/4) / 2 +e4 = (1260 * m_epsilon) ** (1/6) / 2 + +test :: IO () +test = do + print ("e0",e0) + print $ f e0 == 1 + print $ f2 e0 == 1 + -- + print ("e2",e2) + print $ f e2 == f2 e2 + print $ f2 e2 == f4 e2 + -- + print ("e4",e4) + print $ f e4 == f4 e4 + print $ f4 e4 == f6 e4 + diff --git a/math-functions.cabal b/math-functions.cabal --- a/math-functions.cabal +++ b/math-functions.cabal @@ -1,5 +1,5 @@ name: math-functions -version: 0.1.5.2 +version: 0.1.6.0 cabal-version: >= 1.8 license: BSD3 license-file: LICENSE @@ -17,12 +17,13 @@ useful in statistical and numerical computing. extra-source-files: - ChangeLog + changelog.md README.markdown benchmark/*.hs tests/*.hs tests/Tests/*.hs tests/Tests/SpecFunctions/gen.py + doc/sinc.hs library ghc-options: -Wall @@ -38,10 +39,15 @@ Numeric.SpecFunctions Numeric.SpecFunctions.Extra Numeric.Sum + other-modules: + Numeric.SpecFunctions.Internal test-suite tests type: exitcode-stdio-1.0 ghc-options: -Wall -threaded + if arch(i386) + -- The Sum tests require SSE2 on i686 to pass (because of excess precision) + ghc-options: -msse2 hs-source-dirs: tests main-is: tests.hs other-modules: diff --git a/tests/Tests/SpecFunctions.hs b/tests/Tests/SpecFunctions.hs --- a/tests/Tests/SpecFunctions.hs +++ b/tests/Tests/SpecFunctions.hs @@ -78,7 +78,7 @@ ] , 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]] + | n <- [0..1000], k <- [0..n]] ---------------------------------------------------------------- -- Self tests , testProperty "Self-test: 0 <= range01 <= 1" $ \x -> let f = range01 x in f <= 1 && f >= 0 diff --git a/tests/Tests/SpecFunctions_flymake.hs b/tests/Tests/SpecFunctions_flymake.hs new file mode 100644 --- /dev/null +++ b/tests/Tests/SpecFunctions_flymake.hs @@ -0,0 +1,206 @@ +{-# 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 diff --git a/tests/view.hs b/tests/view.hs new file mode 100644 --- /dev/null +++ b/tests/view.hs @@ -0,0 +1,102 @@ +{-# LANGUAGE OverloadedStrings #-} +import Control.Applicative +import Control.Monad +import Numeric.SpecFunctions +import Numeric.MathFunctions.Constants +import CPython.Sugar +import CPython.MPMath +import qualified CPython as Py + +import HEP.ROOT.Plot + + +---------------------------------------------------------------- + + +viewBetaDelta = runPy $ do + addToPythonPath "." + m <- loadMPMath + mpmSetDps m 100 + xs <- forM pqBeta $ \(p,q) -> do x <- fromMPNum =<< mpmLog m =<< mpmBeta m (MPDouble p) (MPDouble q) + return (p,q, relErr x (logBeta p q)) + draws $ do + -- let xs = [ (p,q, logBeta p q `relErr` (logGammaL p + logGammaL q - logGammaL (q+p))) + -- | (p,q) <- pqBeta + -- ] + add $ Graph2D xs + + +pqBeta = [ (p,q) + | p <- logRange 50 0.3 0.6 + , q <- logRange 50 5 6 + ] + where + + + + +viewIBeta x = runPy $ do + addToPythonPath "." + m <- loadMPMath + mpmSetDps m 30 + -- + let n = 40 + let pq = (,) + <$> logRange n 100 1000 + <*> logRange n 100 1000 + -- + xs <- forM pq $ \(p,q) -> do + i <- fromMPNum =<< mpmIncompleteBeta m (MPDouble p) (MPDouble q) (MPDouble x) + return (p,q, incompleteBeta p q x `relErr` i) + -- + draws $ do + add $ Graph2D xs + + +go = runPy $ do + addToPythonPath "." + m <- loadMPMath + mpmSetDps m 16 + -- + print =<< fromMPNum =<< mpmIncompleteBeta m (MPDouble 10) (MPDouble 10) (MPDouble 0.4) + print $ incompleteBeta 10 10 0.4 + + + + +viewLancrox = runPy $ do + addToPythonPath "." + m <- loadMPMath + mpmSetDps m 50 + -- + let xs = logRange 10000 (1e-8) (1e-1) + pl <- forM xs $ \x -> do y0 <- fromMPNum =<< mpmLog m =<< mpmGamma m (MPDouble x) + return (x, y0) + draws $ do + add $ Graph $ [ (x, abs $ y `relErr` logGammaL x) | (x,y) <- pl ] + set $ lineColor RED + -- + add $ Graph $ [ (x, abs $ y `relErr` logGamma x) | (x,y) <- pl ] + set $ lineColor BLUE + -- + set $ xaxis $ logScale ON + -- set $ yaxis $ logScale ON + -- + add $ HLine m_epsilon + add $ HLine $ negate m_epsilon + + +---------------------------------------------------------------- + +relErr :: Double -> Double -> Double +relErr 0 0 = 0 +relErr x y = (x - y) / max (abs x) (abs y) + + + +logRange :: Int -> Double -> Double -> [Double] +logRange n a b + = [ a * r^i | i <- [0 .. n] ] + where + r = (b / a) ** (1 / fromIntegral n) +