## Copyright (C) 2021 Nicholas R. Jankowski ## Copyright (C) 2012 Rik Wehbring ## Copyright (C) 1995-2016 Kurt Hornik ## ## 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; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {} {@var{z} =} binopdf (@var{x}, @var{n}, @var{p}) ## For each element of @var{x}, compute the probability density function ## (PDF) at @var{x} of the binomial distribution with parameters @var{n} ## and @var{p}, where @var{n} is the number of trials and @var{p} is the ## probability of success. ## ## The inputs @var{x}, @var{n}, and @var{p} can be scalars or arrays. If ## more than one is an array, they must be of equal size. The output ## @var{z} will have the same size as the array inputs. ## ## Matlab incompatibility: Octave's @code{binopdf} returns NaN for complex ## input values. ## @end deftypefn ## Author: Nicholas R. Jankowski (Loader version) ## Author: KH (original) ## Description: PDF of the binomial distribution ## Adapted from Loader algorithm for increased calcualtion accuracy. See: ## https://www.r-project.org/doc/reports/CLoader-dbinom-2002.pdf ## and ## http://svn.r-project.org/R/trunk/src/nmath/dbinom.c ## ## p(x,n,p) = p(x,n,x/n)*exp(-D(x,n,p)) ## D(x,n,p) being a deviance term. ## = x ln (x/(np)) + (n-x) ln((n-x)/(n(1-p))) ## ## p(x,n,x/n) rarranges to -> ## sqrt(n/(2pi*x*(n-x))) * exp (d(n)-d(x)-d(n-x)) ## where d(n) = ln(n!*exp(n)/ (sqrt(2pi*n)*n^n)) (overflows) ## = n+gammaln(n+1)-(n+1/2)*log(n)-log(2*pi)/2 (exact) ## = 1./(12 * n) - 1./(360 * n.^3) + 1./(1260 * n.^5) - ## 1./(1680 * n.^7) + 1./(1188 * n.^9);# + O(n^-11); (approx.) function retval = binopdf (x, n, p) if (nargin != 3) print_usage (); endif ## ensure all equal size arrays or scalars, expand scalars to common size. [size_mismatch, x, n, p] = common_size (x, n, p); if (size_mismatch > 0) error ("binopdf: X, N, and P must be of common size or scalars"); endif sz_x = size (x); ## save original size for reshape later x = x(:); n = n(:); p = p(:); ## columns for easier vectorization ## initialize output, preserve class of output if any are singles if (isa (x, "single") || isa (n, "single") || isa (p, "single")); retval = zeros (numel (x), 1, "single"); else retval = zeros (numel (x), 1); endif ## k - index of array locations needing calculation k = (x == fix (x)) & (n == fix (n)) & (n >= 0) & (p >= 0) & (p <= 1) ... & (x >= 0) & (x <= n); nx = n - x; q = 1 - p; ## Catch special cases ahead of calculations: ## Matlab incompatibility: Matalb 2021b returns values for complex inputs ## despite documentation indicating integer and real value inputs required. ## Octave chooses to return an NaN instead. catch_special = (iscomplex (x) | iscomplex (n) | iscomplex (p)); k(catch_special) = false; retval(catch_special) = NaN; ## x = 0 and x = n cases where p != 0 or 1, respectivly ## remove them from k, use alternate calculation to avoid /0 catch_special = (x == 0)& (! catch_special); k(catch_special) = false; retval(catch_special) = exp (n(catch_special) .* log (q(catch_special))); catch_special = (nx == 0) & (! catch_special); k(catch_special) = false; retval(catch_special) = exp (n(catch_special) .* log (p(catch_special))); ## perform Loader pdf calculation on non-trivial elements if (any (k)) retval(k) = loader_expansion (x(k), n(k), p(k), nx(k), q(k)); endif ## Trivial case special outputs: ksp = ((p == 0) & (x == 0)) | (p == 1) & (x == n); retval(ksp) = 1; ## input NaN, n not pos int, or p outside [0,1], ## set output to NaN (overrides 0 or 1) ksp = (n ~= fix (n)) | (n < 0) | (p < 0) | (p > 1) | isnan (x) ... | isnan (n) | isnan (p); retval(ksp) = NaN; retval = reshape (retval, sz_x); ## restore output to input shape endfunction function retval = loader_expansion (x, n, p, nx, q) ## precalculated constants, d_n from n = 0 to 30 ## extended from Loader using octave symbolic vpa ## out to n = 30 d_n = [ 0.08106146679532725821967026359438236013860, 0.04134069595540929409382208140711750802535, 0.02767792568499833914878929274624466659538, 0.02079067210376509311152277176784865633309, 0.01664469118982119216319486537359339114739, 0.01387612882307074799874572702376290856175, 0.01189670994589177009505572411765943862013, 0.01041126526197209649747856713253462919952, 0.00925546218271273291772863663310013611743, 0.00833056343336287125646931865962855220929, 0.00757367548795184079497202421159508389293, 0.00694284010720952986566415266347536265992, 0.00640899418800420706843963108297831257520, 0.00595137011275884773562441604646945832642, 0.00555473355196280137103868995979228464907, 0.00520765591960964044071799685790189865099, 0.00490139594843473786071681819096755442865, 0.00462915374933402859242721316419232323878, 0.00438556024923232426828773634861946570116, 0.00416631969199692245746292338221831613633, 0.00396795421864085961728763680734281467287, 0.00378761806844443457786667706893349200129, 0.00362296022468309470738119836390285473489, 0.00347202138297876696294511542270952959204, 0.00333315563672809287580701911737271025035, 0.00320497022805503801118415655381541759643, 0.00308627868260877706325624133564397946129, 0.00297606398355040882602116255686080370692, 0.00287344936235246638755235148906672207372, 0.00277767492975269360359490376220667282839 ]; stored_dn = numel(d_n); ## indices for precalculated vs to-be-calculated values n_precalc = (n > 0) & (n < stored_dn); x_precalc = (x > 0) & (x < stored_dn); nx_precalc = (nx > 0) & (nx < stored_dn); [delta_n, delta_x, delta_nx] = deal (zeros (size (x))); ## fetch precalculated values delta_n(n_precalc) = d_n(n(n_precalc)); delta_x(x_precalc) = d_n(x(x_precalc)); delta_nx(nx_precalc) = d_n(nx(nx_precalc)); ## calculate any other d(n) values delta_n(!n_precalc) = delta_fn (n(!n_precalc)); delta_x(!x_precalc) = delta_fn (x(!x_precalc)); delta_nx(!nx_precalc) = delta_fn (nx(!nx_precalc)); ## calculate exp(log(pdf)); retval = exp ((delta_n - delta_x - delta_nx + -... deviance (x, n .* p) - ... deviance (nx, n .* q)) - ... 0.5 * (log(2*pi) + log (x) + log (1-x./n))); endfunction function retval = delta_fn (n) ## Stirling formula error term approximations based on Loader paper. ## exact expression, n^n overflows to Inf for n > ~145: ## = log (n!*exp(n)/(sqrt(2pi*n)*n^n)); ## ## Rewritten to avoid overflow out to n> 1e305. accurate to ~10^-12 ## = n + gammaln (n+1) - (n+0.5) * log(n) - log(2*pi)/2; ## ## Approximated as: ## accurate to ~10^-16 for n=30. underflow to 0 at n~10^309 ## = 1/(12n)-1/(360n^3)+1/(1260n^5)- 1/(1680n.^7)+1/(1188n^9) + O(n^-11); ## ## Factored to reduced operation count. Used by Loader and in R: ## 25% faster than unfactored form. ## =(1/12-(1/360-(1/1260-(1/1680-(1/1188)/n^2)/n^2)/n^2)/n^2)/n; nn = n.^2; retval = (0.08333333333333333333333333333333333333333 - ... (0.00277777777777777777777777777777777777778 - ... (0.00079365079365079365079365079365079365079 - ... (0.00059523809523809523809523809523809523810 - ... (0.00084175084175084175084175084175084175084)./nn)./nn)./nn)./nn)./n; endfunction function D = deviance (x, np) ## requires equal length column inputs epsilon = x ./ np; v = (epsilon - 1) ./ (epsilon + 1); vtest = abs (v) < 0.1; if (any (vtest)) ## for abs(v) < 0.1, do taylor expansion for higher precision. Expansion ## term: v^(2j+1)/(2j+1). For abs(v)< 0.1, term drops slowest for max ## abs(v) = 0.1. (n+1)th term is <= 10. jmax = 12; two_jpone = 2 * [1:jmax] + 1; ##sum term 2*j+1 (row vector expansion) D = zeros (numel (epsilon), 1); # D = (x-np)*v + 2*x*sum_over_j(v^2j+1 / 2j+1) D(vtest) = (x(vtest) - np(vtest)) .* v(vtest) + 2 .* x(vtest) .* ... sum (v(vtest).^(two_jpone) ./ two_jpone, 2); D(! vtest) = x(! vtest) .* (log (epsilon(! vtest)) - 1) + np(! vtest); else D = x.* (log (epsilon) - 1) + np; endif endfunction %!shared x,y %! x = [-1 0 1 2 3]; %! y = [0 1/4 1/2 1/4 0]; %!assert (binopdf (x, 2*ones (1,5), 0.5*ones (1,5)), y, eps) %!assert (binopdf (x, 2, 0.5*ones (1,5)), y, eps) %!assert (binopdf (x, 2*ones (1,5), 0.5), y, eps) %!assert (binopdf (x, 2*[0 -1 NaN 1.1 1], 0.5), [0 NaN NaN NaN 0]) %!assert (binopdf (x, 2, 0.5*[0 -1 NaN 3 1]), [0 NaN NaN NaN 0]) %!assert (binopdf ([x, NaN], 2, 0.5), [y, NaN], eps) %!assert (binopdf (cat(3,x,x), 2, 0.5), cat(3,y,y), eps) ## Test Special input values %!assert (binopdf (1, 1, 1), 1) %!assert (binopdf (0, 3, 0), 1) %!assert (binopdf (2, 2, 1), 1) %!assert (binopdf (1, 2, 1), 0) %!assert (binopdf (0, 1.1, 0), NaN) %!assert (binopdf (1, 2, -1), NaN) %!assert (binopdf (1, 2, 1.5), NaN) %!assert (binopdf (i, 2, 0.5), NaN) ##matlab incompatibility %!assert (binopdf (1, i, 0.5), NaN) ##matlab incompatibility %!assert (binopdf (1, 2, i), NaN) ##matlab incompatibility ## Test empty inputs %!assert (binopdf ([], 1, 1), []) %!assert (binopdf (1, [], 1), []) %!assert (binopdf (1, 1, []), []) %!assert (binopdf (ones (1, 0), 2, .5), ones(1, 0)) %!assert (binopdf (ones (0, 1), 2, .5), ones(0, 1)) %!assert (binopdf (ones (0, 1, 2), 2, .5), ones(0, 1, 2)) %!assert (binopdf (1, ones (0, 1, 2), .5), ones(0, 1, 2)) %!assert (binopdf (1, 2, ones (0, 1, 2)), ones(0, 1, 2)) %!assert (binopdf (ones (1, 0, 2), 2, .5), ones(1, 0, 2)) %!assert (binopdf (ones (1, 2, 0), 2, .5), ones(1, 2, 0)) %!assert (binopdf (ones (0, 1, 2), NaN, .5), ones(0, 1, 2)) %!assert (binopdf (ones (0, 1, 2), 2, NaN), ones(0, 1, 2)) ## Test class of input preserved %!assert (binopdf (single ([x, NaN]), 2, 0.5), single ([y, NaN])) %!assert (binopdf ([x, NaN], single (2), 0.5), single ([y, NaN])) %!assert (binopdf ([x, NaN], 2, single (0.5)), single ([y, NaN])) ## Test input validation %!error binopdf () %!error binopdf (1) %!error binopdf (1,2) %!error binopdf (1,2,3,4) %!error binopdf (1, ones (2), ones (3)) %!error binopdf (ones (3), 1, ones (2)) %!error binopdf (ones (3), ones (2), 1) %!error binopdf (ones (3), ones (2), ones (2)) %!error binopdf (ones (2), ones (3), ones (2)) %!error binopdf (ones (2), ones (2), ones (3))