# yiannist/math-functions

Switch branches/tags
Nothing to show
Fetching contributors…
Cannot retrieve contributors at this time
680 lines (634 sloc) 27.3 KB
 -module(spec_functions). -export([ %% * Gamma function log_gamma/1 , log_gamma_l/1 , incomplete_gamma/2 , inv_incomplete_gamma/2 %% * Beta function , log_beta/2 , incomplete_beta/3 , incomplete_beta_/4 , inv_incomplete_beta/3 %% * Logarithm , log1p/1 , log2/1 %% * Factorial , factorial/1 , log_factorial/1 , stirling_error/1 %% * Combinatorics , choose/2 %% * Extras , bd0/2 ]). -export_type([ float64/0 ]). -include("constants.hrl"). %% Special functions and factorials. -type float64() :: float() | 'pos_inf' | 'neg_inf'. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Gamma function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% %% Adapted from http://people.sc.fsu.edu/~burkardt/f_src/asa245/asa245.html %% %% Compute the logarithm of the gamma function Gamma(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 'pos_inf' if the input is outside of the range %% (0 < x <= 1e305). -spec log_gamma(float()) -> float64(). log_gamma(X) -> %% A bunch of constants: Y = math:log(X), Alr2pi = 0.918938533204673, K = X * (Y - 1) - 0.5 * Y + Alr2pi, 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, %% Use a helper instead of guards! ;-) LogGamma_wrap = fun (X00) -> {A, B, C} = case X00 < 0.5 of true -> {-Y, X00 + 1, X00}; false -> {0, X00, X00 - 1} end, LogGamma = fun (X0) when X0 =< 0 -> 'pos_inf'; (X0) when X0 < 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); (X0) when X0 < 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); (X0) when X0 < 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); (X0) when X0 > 5.1e5 -> K; (_) -> K + X1 * ((R4_2 * X2 + R4_1) * X2 + R4_0) / ((X2 + R4_4) * X2 + R4_3) end, LogGamma(X00) end, LogGamma_wrap(X). %% Compute the logarithm of the gamma function, G(x). Uses a %% Lanczos approximation. %% %% This function is slower than 'log_gamma', but gives 14 or more %% significant decimal digits of accuracy, except around x = 1 and %% x = 2, where the function goes to zero. %% %% Returns M_HUGE; if the input is outside of the range %% (0 < x <= 1e305). -spec log_gamma_l(float()) -> float64(). log_gamma_l(X) when X =< 0 -> 'pos_inf'; log_gamma_l(X) -> A0 = 0.9999999999995183, X65 = X + 6.5, Go = fun (K, {L, T}) -> {L + K / T, T - 1} end, A = [ 0.1659470187408462e-06, 0.9934937113930748e-05 , -0.1385710331296526, 12.50734324009056 , -176.6150291498386, 771.3234287757674 , -1259.139216722289, 676.5203681218835 ], Fini = fun ({L, _}) -> math:log(L + A0) + math:log(?M_SQRT_2_PI) - X65 + (X - 0.5) * math:log(X65) end, Fini(lists:foldl(Go, {0, X + 7}, A)). %% Compute the log gamma correction factor for X > 10. This %% correction factor is suitable for an alternate (but less %% numerically accurate) definition of 'log_gamma': %% %% > lgg x = 0.5 * log(2*pi) + (x-0.5) * log x - x + logGammaCorrection x -spec log_gamma_correction(float()) -> float(). log_gamma_correction(X) when X < 10 -> exit({?MODULE, log_gamma_correction, "NaN"}); log_gamma_correction(X) -> Big = 94906265.62425156, T = 10 / X, Coeffs = [ 0.1666389480451863247205729650822e+0 , -0.1384948176067563840732986059135e-4 , 0.9810825646924729426157171547487e-8 , -0.1809129475572494194263306266719e-10 , 0.6221098041892605227126015543416e-13 , -0.3399615005417721944303330599666e-15 , 0.2683181998482698748957538846666e-17 ], case X < Big of true -> polynomial:chebyshev_broucke(T * T * 2 - 1, Coeffs) / X; false -> 1 / (X * 12) end. %% Compute the normalized lower incomplete gamma function %% g(s, x). Normalization means that g(s, inf)=1. %% Uses Algorithm AS 239 by Shea. -spec incomplete_gamma(float(), float()) -> float64(). incomplete_gamma(P0, X0) -> Norm = fun (A) -> 0.5 * math:erfc( -A / ?M_SQRT_2 ) end, Limit = -88, Tolerance = 1.0e-14, Overflow = 1.0e37, CF = fun () -> A = 1 - P0, B = A + X0 + 1, P3 = X0 + 1, P4 = X0 * B, cont_frac(Tolerance, Overflow, A, B, 0, 1.0, X0, P3, P4, P3 / P4) end, IncompleteGamma = fun (P, X) when ((X < 0) or (P =< 0)) -> 'pos_inf'; (_, 0) -> 0; (P, X) when P >= 1000 -> Norm(3 * math:sqrt(P) * (math:pow(X / P, 1 / 3) + 1 / (9 * P) - 1)); (_, X) when X > 1.0e8 -> 1; (P, X) when ((X =< 1) or (X < P)) -> A = P * math:log(X) - X - log_gamma(P + 1), G = A + math:log(pearson(Tolerance, X, P, 1.0, 1.0)), case G > Limit of true -> math:exp(G); false -> 0 end; (P, X) -> G = P*math:log(X) - X - log_gamma(P) + math:log(CF()), case G > Limit of true -> 1 - math:exp(G); false -> 1 end end, IncompleteGamma(P0, X0). -spec pearson(float(), float(), float(), float(), float()) -> float(). pearson(Tolerance, X, A, C, G) -> A2 = A + 1, C2 = C * X / A2, G2 = G + C2, case C2 =< Tolerance of true -> G2; false -> pearson(Tolerance, X, A2, C2, G2) end. -spec cont_frac(float(), float(), float(), float(), integer(), float() , float(), float(), float(), float()) -> float(). cont_frac(Tolerance, Overflow, A, B, C, P1, P2, P3, P4, G) -> A2 = A + 1, B2 = B + 2, C2 = C + 1, AN = A2 * C2, P5 = B2 * P3 - AN * P1, P6 = B2 * P4 - AN * P2, RN = P5 / P6, F = fun (N) -> case abs(P5) > Overflow of true -> N / Overflow; false -> N end end, case abs(G - RN) =< min(Tolerance, Tolerance * RN) of true -> G; false -> cont_frac(Tolerance, Overflow, A2, B2, C2, F(P3), F(P4), F(P5), F(P6), RN) end. %% Adapted from Numerical Recipes Â§6.2.1 %% Inverse incomplete gamma function. It's approximately inverse of %% 'incomplete_gamma' for the same s. So following equality %% approximately holds: %% %% > invIncompleteGamma s . incompleteGamma s = id %% %% For inv_incomplete_gamma s p, s must be positive and p must be %% in [0,1] range. %%STUB: Might notice some deviation from Haskell results (~1e-11). -spec inv_incomplete_gamma(float(), float()) -> float64(). inv_incomplete_gamma(A, P) -> %% Calculate initial guess for root Guess = fun () when A > 1 -> Lg = case P < 0.5 of true -> P; false -> 1 - P end, T = math:sqrt(-2 * math:log(Lg)), X1 = (2.30753 + T * 0.27061) / (1 + T * (0.99229 + T * 0.04481)) - T, X2 = case P < 0.5 of true -> -X1; false -> X1 end, max(1.0e-3, (A * math:pow((1 - 1/(9 * A) - X2 / (3 * math: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 () -> T = 1 - A * (0.253 + A * 0.12), case P < T of true -> math:pow(P / T, 1 / A); false -> 1 - math:log( 1 - (P - T) / (1 - T) ) end end, Fun = fun () when A =< 0 -> exit({?MODULE, inv_incomplete_gamma, "'a' must be positive!"}); () when ((P < 0) or (P > 1)) -> exit({?MODULE, inv_incomplete_gamma, "'p' must be in [0,1]"}); () when P == 0 -> 0; () when P == 1 -> 'pos_inf'; () -> inv_incomplete_gamma__go({A, P}, 0, Guess()) end, Fun(). %% Solve equation γ(a,x) = p using Halley method -spec inv_incomplete_gamma__go({float(), float()}, integer(), float()) -> float(). inv_incomplete_gamma__go(_, I, X) when I >= 12 -> X; inv_incomplete_gamma__go({A, P}, I, X) -> %% Value of γ(a,x) - p F = incomplete_gamma(A, X) - P, %% Constants EPS = 1.0e-8, A1 = A - 1, GLN = log_gamma(A), %% dγ(a,x)/dx F2 = case A > 1 of true -> LnA1 = math:log(A1), AFAC = math:exp(A1 * (LnA1 - 1) - GLN), AFAC * math:exp( -(X - A1) + A1 * (math:log(X) - LnA1)); false -> math:exp( -X + A1 * math:log(X) - GLN) end, U = F / F2, %% Halley correction to Newton-Rapson step Corr = U * (A1 / X - 1), DX = U / (1 - 0.5 * min(1.0, Corr)), %% New approximation to x X2 = case X < DX of true -> 0.5 * X; %% Do not go below 0 false -> X - DX end, case abs(DX) < EPS * X2 of true -> X2; false -> inv_incomplete_gamma__go({A, P}, I + 1, X2) end. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Beta function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Compute the natural logarithm of the beta function. -spec log_beta(float(), float()) -> float64(). log_beta(A, B) -> P = min(A, B), Q = max(A, B), F = fun () when P < 0 -> exit({?MODULE, log_beta, "NaN"}); () when P == 0 -> 'pos_inf'; () when P >= 10 -> PQ = P + Q, PPQ = P / PQ, C = log_gamma_correction(Q) - log_gamma_correction(PQ), math:log(Q) * (-0.5) + ?M_LN_SQRT_2_PI + log_gamma_correction(P) + C + (P-0.5) * math:log(PPQ) + Q * log1p(-PPQ); () when Q >= 10 -> PQ = P + Q, PPQ = P / PQ, C = log_gamma_correction(Q) - log_gamma_correction(PQ), log_gamma(P) + C + P - P * math:log(PQ) + (Q - 0.5) * log1p(-PPQ); () -> log_gamma(P) + log_gamma(Q) - log_gamma(P + Q) end, F(). %% Regularized incomplete beta function. Uses algorithm AS63 by %% Majumder and Bhattachrjee. -spec incomplete_beta(float(), float(), float()) -> float(). incomplete_beta(P, Q, X) -> incomplete_beta_(log_beta(P, Q), P, Q, X). %% Regularized incomplete beta function. Same as 'incomplete_beta' %% but also takes logarithm of beta function as parameter. -spec incomplete_beta_(float(), float(), float(), float()) -> float(). incomplete_beta_(_, P, Q, _) when ((P =< 0) or (Q =< 0)) -> exit({?MODULE, incomplete_beta_, "p <= 0 || q <= 0"}); incomplete_beta_(_, _, _, X) when ((X < 0) or (X > 1)) -> exit({?MODULE, incomplete_beta_, "x out of [0,1] range"}); incomplete_beta_(_, _, _, X) when ((X == 0) or (X == 1)) -> X; incomplete_beta_(Beta, P, Q, X) when P >= (P + Q) * X -> incomplete_beta_worker(Beta, P, Q, X); incomplete_beta_(Beta, P, Q, X) -> 1 - incomplete_beta_worker(Beta, Q, P, 1 - X). %% Worker for incomplete beta function. It is separate function to %% avoid confusion with parameter during parameter swapping -spec incomplete_beta_worker(float(), float(), float(), float()) -> float(). incomplete_beta_worker(Beta, P, Q, X) -> incomplete_beta_worker__go({Beta, P, Q, X}, P + Q, trunc(Q + (1 - X) * (P + Q)), 1, 1.0, 1.0). %% Loop. -spec incomplete_beta_worker__go({float(), float(), float(), float()}, float(), integer(), integer(), float(), float()) -> float(). incomplete_beta_worker__go({Beta, P, Q, X}, PSQ, NS, AI, Term, Betain) -> %% New values EPS = 1.0e-15, CX = 1 - X, Fact = case NS > 0 of true -> (Q - AI) * X / CX; false -> case NS == 0 of true -> (Q - AI) * X; false -> PSQ * X end end, Term2 = Term * Fact / (P + AI), Betain2 = Betain + Term2, %% Iterations are complete DB = abs(Term2), Done = (DB =< EPS) and (DB =< EPS * Betain2), PSQ2 = case NS < 0 of true -> PSQ + 1; false -> PSQ end, case Done of true -> Betain2 * math:exp(P*math:log(X) + (Q-1) * math:log(CX) - Beta) / P; false -> incomplete_beta_worker__go({Beta, P, Q, X}, PSQ2, NS-1 , AI+1, Term2, Betain2) end. %% Compute inverse of regularized incomplete beta function. Uses %% initial approximation from AS109 and Halley method to solve equation. -spec inv_incomplete_beta(float(), float(), float()) -> float(). inv_incomplete_beta(P, Q, _) when ((P =< 0) or (Q =< 0)) -> exit({?MODULE, inv_incomplete_beta, "p <=0 || q <=0"}); inv_incomplete_beta(_, _, A) when ((A < 0) or (A > 1)) -> exit({?MODULE, inv_incomplete_beta, "bad a"}); inv_incomplete_beta(_, _, A) when ((A == 0) or (A == 1)) -> A; inv_incomplete_beta(P, Q, A) when A > 0.5 -> 1 - inv_incomplete_beta_worker(log_beta(P, Q), Q, P, 1 - A); inv_incomplete_beta(P, Q, A) -> inv_incomplete_beta_worker(log_beta(P, Q), P, Q, A). -spec inv_incomplete_beta_worker(float(), float(), float(), float()) -> float(). inv_incomplete_beta_worker(Beta, P, Q, A) -> %% Calculate initial guess R = math:sqrt( -math:log(A * A) ), Y = R - ( 2.30753 + 0.27061 * R ) / ( 1.0 + ( 0.99229 + 0.04481 * R ) * R ), T = 1 / (9 * Q), T2 = 2 * Q * math:pow(1 - T + Y * math:sqrt(T), 3), T3 = (4*P + 2*Q - 2) / T2, Guess = fun () when ((P > 1) and (Q > 1)) -> RR = (Y*Y - 3) / 6, SS = 1 / (2*P - 1), TT = 1 / (2*Q - 1), HH = 2 / (SS + TT), WW = Y * math:sqrt(HH + RR) / HH - (HH - SS) * (RR + 5/6 - 2 / (3 * HH)), P / (P + Q * math:exp(2 * WW)); () when T2 =< 0 -> 1 - math:exp( (math:log((1 - A) * Q) + Beta) / Q ); () when T3 =< 1 -> math:exp( (math:log(A * P) + Beta) / P ); () -> 1 - 2 / (T3 + 1) end, inv_incomplete_beta_worker__go({Beta, P, Q, A}, 0, Guess()). -spec inv_incomplete_beta_worker__go({float(), float(), float(), float()}, integer(), float()) -> float(). inv_incomplete_beta_worker__go({Beta, P, Q, A}, I, X) -> Fun = fun () when ((X == 0) or (X == 1) or (I >= 10)) -> X; () -> P1 = P - 1, Q1 = Q - 1, F = incomplete_beta_(Beta, P, Q, X) - A, F2 = math:exp(P1 * math:log(X) + Q1 * math:log(1-X) - Beta), U = F / F2, DX = U / (1 - 0.5 * min(1, (U * (P1 / X - Q1 / (1 - X))))), Z = X - DX, X2 = fun () when Z < 0 -> X / 2; () when Z > 1 -> (X + 1) / 2; () -> Z end, case abs(DX) =< 16 * ?M_EPSILON * X of true -> X; false -> inv_incomplete_beta_worker__go({Beta, P, Q, A} , I + 1, X2()) end end, Fun(). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 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. -spec log1p(float()) -> float64(). log1p(0) -> 0; log1p(-1) -> 'neg_inf'; log1p(X) when X < -1 -> exit({?MODULE, log1p, "NaN"}); log1p(X) -> X2 = abs(X), Coeffs = [ 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 ], Fun = fun () when X2 < ?M_EPSILON * 0.5 -> X; () when ((X >= 0) and (X < 1.0e-8)) -> X * (1 - X * 0.5); () when X2 < 0.375 -> X * (1 - X * polynomial:chebyshev_broucke(X / 0.375, Coeffs)); () -> math:log(1 + X) end, Fun(). %% Compute the logarithm in base 2 of the given value. -spec log2(integer()) -> integer(). log2(V0) when V0 =< 0 -> exit({?MODULE, log2, "invalid input"}); log2(V0) -> log2__go(5, 0, V0). -spec log2__go(integer(), integer(), integer()) -> integer(). log2__go(-1, R, _) -> R; log2__go( I, R, V) -> BV = [16#2, 16#c, 16#f0, 16#ff00, 16#ffff0000, 16#ffffffff00000000], SV = [1,2,4,8,16,32], case (V band lists:nth(I+1, BV)) =/= 0 of true -> SI = lists:nth(I+1, SV), log2__go(I-1, R bor SI, V bsr SI); false -> log2__go(I-1, R, V) end. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Factorial %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Compute the factorial function n!. Returns ?M_HUGE if the %% input is above 170 (above which the result cannot be represented by %% a 64-bit 'float'). -spec factorial(integer()) -> float64(). factorial(N) when N < 0 -> exit({?MODULE, factorial, "negative input"}); factorial(N) when N =< 1 -> 1.0; factorial(N) when N =< 170 -> Prod = fun (X, Acc) -> X * Acc end, lists:foldl(Prod, 1.0, lists:seq(2, N)); factorial(_) -> 'pos_inf'. %% Compute the natural logarithm of the factorial function. Gives %% 16 decimal digits of precision. -spec log_factorial(integer()) -> float(). log_factorial(N) when N =< 14 -> math:log(factorial(N)); log_factorial(N) -> X = N + 1, Y = 1 / (X * X), Z = ((-(5.95238095238e-4 * Y) + 7.936500793651e-4) * Y - 2.7777777777778e-3) * Y + 8.3333333333333e-2, (X - 0.5) * math:log(X) - X + 9.1893853320467e-1 + Z / X. %% 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) -spec stirling_error(float()) -> float(). stirling_error(N) -> NN = N * N, 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 = [ 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 ], Fun = fun () when N =< 15 -> ProperFraction = fun (N0) -> (2 * N0) == (2 * trunc(N0)) end, case ProperFraction(N) of true -> lists:nth(2 * trunc(N), Sfe); false -> log_gamma(N+1.0) - (N+0.5) * math:log(N) + N - ?M_LN_SQRT_2_PI end; () when N > 500 -> (S0-S1/NN)/N; () when N > 80 -> (S0-(S1-S2/NN)/NN)/N; () when N > 35 -> (S0-(S1-(S2-S3/NN)/NN)/NN)/N; () -> (S0-(S1-(S2-(S3-S4/NN)/NN)/NN)/NN)/N end, Fun(). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Combinatorics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Quickly compute the natural logarithm of n `choose` k, with %% no checking. -spec log_choose_fast(float(), number()) -> float(). log_choose_fast(N, K) -> -math:log(N + 1) - log_beta(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: %% %% 1> choose(7, 3). %% 35 -spec choose(integer(), integer()) -> float(). choose(N, K) when K > N -> 0; choose(N, K) -> K2 = min(K, N-K), NK = N - K2, Approx = math:exp(log_choose_fast(N, K2)), Fun = fun () when K2 < 50 -> Go = fun (I, A) -> A * (NK+I) / I end, lists:foldl(Go, 1, lists:seq(1, K2)); () when Approx < ?MAX64 -> round(Approx); %STUB: Is ?MAX64 correct? () -> Approx end, Fun(). %% From Numeric.SpecFunctions.Extra: %% Evaluate the deviance term: x log(x/np) + np - x. %% STUB: stripped isInfinite guards! -spec bd0(float(), float()) -> float(). bd0(_, 0) -> exit({?MODULE, bd0, "NaN"}); bd0(X, NP) -> X_NP = X - NP, case abs(X_NP) >= 0.1 * (X + NP) of true -> X * math:log(X/NP) - X_NP; false -> V = X_NP / (X + NP), S0 = X_NP * V, EJ0 = 2 * X * V, VV = V * V, bd0__go(VV, 1, EJ0 * VV, S0) end. -spec bd0__go(float(), pos_integer(), float(), float()) -> float(). bd0__go(VV, J, EJ, S) -> S2 = S + EJ / (2 * J + 1), case S2 == S of true -> S2; %% FIXME: Comparing floats for equality! false -> bd0__go(VV, J + 1, (EJ * VV), S2) end. %% References: %% %% * Lanczos, C. (1964) A precision approximation of the gamma %% function. /SIAM Journal on Numerical Analysis B/ %% 1:86–96. %% %% * Loader, C. (2000) Fast and Accurate Computation of Binomial %% Probabilities. %% %% * Macleod, A.J. (1989) Algorithm AS 245: A robust and reliable %% algorithm for the logarithm of the gamma function. %% /Journal of the Royal Statistical Society, Series C (Applied Statistics)/ %% 38(2):397–402. %% %% * Shea, B. (1988) Algorithm AS 239: Chi-squared and incomplete %% gamma integral. /Applied Statistics/ %% 37(3):466–473. %% %% * Majumder, K.L., Bhattacharjee, G.P. (1973) Algorithm AS 63: The %% Incomplete Beta Integral. /Journal of the Royal Statistical %% Society. Series C (Applied Statistics)/ Vol. 22, No. 3 (1973), %% pp. 409-411. %% %% * Majumder, K.L., Bhattacharjee, G.P. (1973) Algorithm AS 64: %% Inverse of the Incomplete Beta Function Ratio. /Journal of the %% Royal Statistical Society. Series C (Applied Statistics)/ %% Vol. 22, No. 3 (1973), pp. 411-414 %% %% %% * Cran, G.W., Martin, K.J., Thomas, G.E. (1977) Remark AS R19 %% and Algorithm AS 109: A Remark on Algorithms: AS 63: The %% Incomplete Beta Integral AS 64: Inverse of the Incomplete Beta %% Function Ratio. /Journal of the Royal Statistical Society. Series %% C (Applied Statistics)/ Vol. 26, No. 1 (1977), pp. 111-114 %%