module Numeric.SpecFunctions (
logGamma
, logGammaL
, incompleteGamma
, invIncompleteGamma
, logBeta
, incompleteBeta
, incompleteBeta_
, invIncompleteBeta
, log1p
, log2
, factorial
, logFactorial
, stirlingError
, choose
) where
import Data.Bits ((.&.), (.|.), shiftR)
import Data.Int (Int64)
import Data.Word (Word64)
import Data.Number.Erf (erfc)
import qualified Data.Vector.Unboxed as U
import Numeric.Polynomial.Chebyshev (chebyshevBroucke)
import Numeric.MathFunctions.Constants (m_epsilon, m_sqrt_2_pi, m_ln_sqrt_2_pi,
m_NaN, m_neg_inf, m_pos_inf, m_sqrt_2)
logGamma :: Double -> Double
logGamma x
| x <= 0 = m_pos_inf
| 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 > 5.1e5 = 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 * (y1) 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 !Double !Double
logGammaL :: Double -> Double
logGammaL x
| x <= 0 = m_pos_inf
| otherwise = fini . U.foldl' go (L 0 (x+7)) $ a
where fini (L l _) = log (l+a0) + log m_sqrt_2_pi x65 + (x0.5) * log x65
go (L l t) k = L (l + k / t) (t1)
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
]
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
]
incompleteGamma :: Double
-> Double
-> Double
incompleteGamma p x
| x < 0 || p <= 0 = m_pos_inf
| x == 0 = 0
| p >= 1000 = norm (3 * sqrt p * ((x/p) ** (1/3) + 1/(9*p) 1))
| 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
norm a = erfc ( a / m_sqrt_2)
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
invIncompleteGamma :: Double -> Double -> Double
invIncompleteGamma a p
| a <= 0 =
error $ "Statistics.Math.invIncompleteGamma: a must be positive. Got: " ++ show a
| p < 0 || p > 1 =
error $ "Statistics.Math.invIncompleteGamma: p must be in [0,1] range. Got: " ++ show p
| p == 0 = 0
| p == 1 = 1 / 0
| otherwise = loop 0 guess
where
loop :: Int -> Double -> Double
loop i x
| i >= 12 = x
| otherwise =
let
f = incompleteGamma a x p
f' | a > 1 = afac * exp( (x a1) + a1 * (log x lna1))
| otherwise = exp( x + a1 * log x gln)
u = f / f'
corr = u * (a1 / x 1)
dx = u / (1 0.5 * min 1.0 corr)
x' | x < dx = 0.5 * x
| otherwise = x dx
in if abs dx < eps * x'
then x'
else loop (i+1) x'
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)
| otherwise =
let t = 1 a * (0.253 + a*0.12)
in if p < t
then (p / t) ** (1 / a)
else 1 log( 1 (pt) / (1t))
a1 = a 1
lna1 = log a1
afac = exp( a1 * (lna1 1) gln )
gln = logGamma a
eps = 1e-8
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
incompleteBeta :: Double
-> Double
-> Double
-> Double
incompleteBeta p q = incompleteBeta_ (logBeta p q) p q
incompleteBeta_ :: Double
-> Double
-> Double
-> Double
-> Double
incompleteBeta_ beta p q x
| p <= 0 || q <= 0 = error "p <= 0 || q <= 0"
| x < 0 || x > 1 = error "x < 0 || x > 1"
| x == 0 || x == 1 = x
| p >= (p+q) * x = incompleteBetaWorker beta p q x
| otherwise = 1 incompleteBetaWorker beta q p (1 x)
incompleteBetaWorker :: Double -> Double -> Double -> Double -> Double
incompleteBetaWorker beta p q x = loop (p+q) (truncate $ q + cx * (p+q) :: Int) 1 1 1
where
eps = 1e-15
cx = 1 x
loop psq ns 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
term' = term * fact / (p + ai)
betain' = betain + term'
fact | ns > 0 = (q ai) * x/cx
| ns == 0 = (q ai) * x
| otherwise = psq * x
done = db <= eps && db <= eps*betain' where db = abs term'
psq' = if ns < 0 then psq + 1 else psq
invIncompleteBeta :: Double
-> Double
-> Double
-> Double
invIncompleteBeta p q a
| p <= 0 || q <= 0 = error "p <= 0 || q <= 0"
| a < 0 || a > 1 = error "bad 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
invIncompleteBetaWorker beta p q a = loop (0::Int) guess
where
p1 = p 1
q1 = q 1
loop !i !x
| x == 0 || x == 1 = x
| i >= 10 = x
| abs dx <= 16 * m_epsilon * x = x
| otherwise = loop (i+1) x'
where
f = incompleteBeta_ beta p q x a
f' = exp $ p1 * log x + q1 * log (1 x) beta
u = f / f'
dx = u / (1 0.5 * min 1 (u * (p1 / x q1 / (1 x))))
x' | z < 0 = x / 2
| z > 1 = (x + 1) / 2
| otherwise = z
where z = x dx
guess
| p > 1 && q > 1 =
let rr = (y*y 3) / 6
ss = 1 / (2*p 1)
tt = 1 / (2*q 1)
hh = 2 / (ss + tt)
ww = y * sqrt(hh + rr) / hh (tt ss) * (rr + 5/6 2 / (3 * hh))
in p / (p + q * exp(2 * ww))
| t' <= 0 = 1 exp( (log((1 a) * q) + beta) / q )
| t'' <= 1 = exp( (log(a * p) + beta) / p )
| otherwise = 1 2 / (t'' + 1)
where
r = sqrt ( log ( a * a ) )
y = r ( 2.30753 + 0.27061 * r )
/ ( 1.0 + ( 0.99229 + 0.04481 * r ) * r )
t = 1 / (9 * q)
t' = 2 * q * (1 t + y * sqrt t) ** 3
t'' = (4*p + 2*q 2) / t'
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
]
log2 :: Int -> Int
log2 v0
| v0 <= 0 = error "Statistics.Math.log2: invalid input"
| 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 (i1) (r .|. si) (v `shiftR` si)
| otherwise = go (i1) 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 :: Int -> Double
factorial n
| n < 0 = error "Statistics.Math.factorial: negative input"
| n <= 1 = 1
| n <= 170 = U.product $ U.map fromIntegral $ U.enumFromTo 2 n
| otherwise = m_pos_inf
logFactorial :: Int -> Double
logFactorial n
| n <= 14 = log (factorial 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
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 = (s0s1/nn)/n
| n > 80 = (s0(s1s2/nn)/nn)/n
| n > 35 = (s0(s1(s2s3/nn)/nn)/nn)/n
| otherwise = (s0(s1(s2(s3s4/nn)/nn)/nn)/nn)/n
where
nn = n*n
s0 = 0.083333333333333333333
s1 = 0.00277777777777777777778
s2 = 0.00079365079365079365079365
s3 = 0.000595238095238095238095238
s4 = 0.0008417508417508417508417508
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 ]
logChooseFast :: Double -> Double -> Double
logChooseFast n k = log (n + 1) logBeta (n k + 1) (k + 1)
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 (nk)
approx = exp $ logChooseFast (fromIntegral n) (fromIntegral k')
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