## Copyright (C) 2016 Marco Caliari, Nir Krakauer, Andrew Lachlan ## ## This program is free software; you can redistribute it and/or modify it under ## the terms of the GNU General Public License as published by the Free Software ## Foundation; either version 3 of the License, or (at your option) any later ## version. ## ## This program is distributed in the hope that it will be useful, but WITHOUT ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more ## details. ## ## You should have received a copy of the GNU General Public License along with ## this program; if not, see . function y = D(a,x) ## Compute exp(-x)*x^a/Gamma(a+1) in a stable way for x and a large ## ## L. Knusel, Computation of the Chi-square and Poisson distribution, ## SIAM J. Sci. Stat. Comput., 7(3), 1986 ## which quotes Section 5, Abramowitz&Stegun 6.1.40, 6.1.41 athresh = 10; ## FIXME: can be better tuned? if (x == 0) y = 0; else if ((x > 0) && (a > athresh)) lnGa = log (2 * pi *a) / 2+1 / (12 * a) - 1 / (360 * a ^ 3) + ... 1 / (1260 * a ^ 5) - 1 / (1680 * a ^ 7) + 1 / (1188 * a ^ 9)- ... 691 / (87360 * a ^ 11) + 1 / (156 * a ^ 13) - ... 3617 / (122400 * a ^ 15) + 43867 / (244188 * a ^ 17) - ... 174611 / (125400 * a ^ 19); if (a >= x) lns = log1p ((a - x) / x); else lns = -log1p ((x - a) / a); endif y = exp ((a - x) - a * lns - lnGa); else ## standard formula for a not so large y = exp(a * log (x) - x - gammaln (a + 1)); endif if ((x < 0) && (a == fix (a))) % remove spurious imaginary part y = real (y); endif endif %!test %! assert(D(10,10),0.12511003572113329898476497894772544708420990097708600,-7*eps) %%!test %%! assert(D(100,10),4.864649182067610436458701223522755117545181424119392e-63,-115*eps) %!test %! assert(D(100,100),0.0398609968091471352339206494591387451322894644246,-3*eps) %!test %! assert(D(200,200),0.0281977276859208217987020623393630120002075706645,-3*eps) %!test %! assert(D(200,201),0.0281275549310563187334442096017465724513046028959,-3*eps) %!test %! assert(D(1000,1000),0.01261461134872149971803693647457875764715489427,-4*eps) %!test %! assert(D(10,-10),6.069903492836947893782490257188118487200593203967e7,-8*eps) %!test %! assert(D(10,-10.1),7.410115412614356278962877535200296949490313010e7,-20*eps) %!test %! assert(D(9,-10),-6.069903492836947893782490257188118487200593203967e7,-8*eps) %!test %! assert(D(10.1,-10),5.741783267595397392122436851864044460198744897912388e7+1i*1.865618475142280818086513704349404376045316209289692e7,-25*eps)