## Copyright (C) 2012 Rik Wehbring
## Copyright (C) 1995-2015 Kurt Hornik
## Copyright (C) 2015 Lachlan Andrew
## Copyright (C) 2014 Mike Giles
##
## This file is part of Octave.
##
## Octave 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.
##
## Octave 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 Octave; see the file COPYING. If not, see
## .
## -*- texinfo -*-
## @deftypefn {Function File} {} poissinv (@var{x}, @var{lambda})
## For each element of @var{x}, compute the quantile (the inverse of the CDF)
## at @var{x} of the Poisson distribution with parameter @var{lambda}.
## Both @var{x} and @var{lambda} may be either scalars or vectors. If
## they are vectors they must have the same size.
## @end deftypefn
## Author: Lachlan
## based on code by
## KH
## Mike Giles
## Description: Quantile function of the Poisson distribution
function inv = poissinv (x, lambda)
if (nargin != 2)
print_usage ();
endif
if (! isscalar (lambda))
[retval, x, lambda] = common_size (x, lambda);
if (retval > 0)
error ("poissinv: X and LAMBDA must be of common size or scalars");
endif
lam = lambda;
else
lam = repmat (lambda, size (x)); # vector lambda, for analytic_approx
endif
if (iscomplex (x) || iscomplex (lambda))
error ("poissinv: X and LAMBDA must not be complex");
endif
if (isa (x, "single") || isa (lambda, "single"))
inv = zeros (size (x), "single");
else
inv = zeros (size (x));
endif
k = (x < 0) | (x > 1) | isnan (x) | !(lambda > 0);
inv(k) = NaN;
k = (x == 1) & (lambda > 0);
inv(k) = Inf;
k = ((x > 0) & (x < 1) & (lambda > 0));
# Find expression Mike Giles's magic when lambda isn't small
kk = k & (lambda >= 4);
if (any (kk))
inv(kk) = analytic_approx (x(kk), lam(kk));
endif
# Use bottom-up search for cases when the answer is small, or not set above
k = k & (inv(k) < 10);
k = find (k);
if (! isempty (k))
if (isscalar (lambda))
cdf = exp (-lambda) * ones (size (k));
else
cdf = exp (-lambda(k));
endif
inv(k) = 0;
for i = 1:100 ## Truncate at 100 for lambda<4
m = find (cdf < x(k));
if (any (m))
inv(k(m)) += 1;
if (isscalar (lambda))
cdf = cdf + poisspdf (i, lambda);
else
cdf(m) = cdf(m) + poisspdf (i, lambda(k(m)));
endif
else
break;
endif
endfor
# if the truncation kicked in, use the analytic approximation even if
# lambda < 4.
if (any (m))
inv(k(m)) = analytic_approx (x(k(m)), lam(k(m)));
endif
endif
endfunction
## The following is based on Mike Giles's CUDA implementation,
## [http://people.maths.ox.ac.uk/gilesm/codes/poissinv/poissinv_cuda.h]
## which is copyright by the University of Oxford
## and is provided under the terms of the GNU GPLv3 license:
## http://www.gnu.org/licenses/gpl.html
function inv = analytic_approx (x, lambda)
s = norminv (x, 0, 1) ./ sqrt (lambda);
k = ((s > -0.6833501) & (s < 1.777993));
# use polynomial approximations in central region
if (any (k))
lam = lambda(k);
if (isscalar (s))
sk = s;
else
sk = s(k);
endif
# polynomial approximation to f^{-1}(s) - 1
rm = 2.82298751e-07;
rm = -2.58136133e-06 + rm.*sk;
rm = 1.02118025e-05 + rm.*sk;
rm = -2.37996199e-05 + rm.*sk;
rm = 4.05347462e-05 + rm.*sk;
rm = -6.63730967e-05 + rm.*sk;
rm = 0.000124762566 + rm.*sk;
rm = -0.000256970731 + rm.*sk;
rm = 0.000558953132 + rm.*sk;
rm = -0.00133129194 + rm.*sk;
rm = 0.00370367937 + rm.*sk;
rm = -0.0138888706 + rm.*sk;
rm = 0.166666667 + rm.*sk;
rm = sk + sk.*(rm.*sk);
# polynomial approximation to correction c0(r)
t = 1.86386867e-05;
t = -0.000207319499 + t.*rm;
t = 0.0009689451 + t.*rm;
t = -0.00247340054 + t.*rm;
t = 0.00379952985 + t.*rm;
t = -0.00386717047 + t.*rm;
t = 0.00346960934 + t.*rm;
t = -0.00414125511 + t.*rm;
t = 0.00586752093 + t.*rm;
t = -0.00838583787 + t.*rm;
t = 0.0132793933 + t.*rm;
t = -0.027775536 + t.*rm;
t = 0.333333333 + t.*rm;
# O(1/lam) correction
y = -0.00014585224;
y = 0.00146121529 + y.*rm;
y = -0.00610328845 + y.*rm;
y = 0.0138117964 + y.*rm;
y = -0.0186988746 + y.*rm;
y = 0.0168155118 + y.*rm;
y = -0.013394797 + y.*rm;
y = 0.0135698573 + y.*rm;
y = -0.0155377333 + y.*rm;
y = 0.0174065334 + y.*rm;
y = -0.0198011178 + y.*rm;
y /= lam;
inv(k) = floor (lam + (y+t)+lam.*rm);
endif
k = !k & s > -sqrt (2);
if (any (k))
# Newton iteration
r = 1 + s(k);
r2 = r + 1;
while (abs (r - r2) > 1e-5)
t = log (r);
r2 = r;
s2 = sqrt (2 * ((1-r) + r.*t));
s2(r<1) *= -1;
r = r2 - (s2 - s(k)) .* s2 ./ t;
if (r < 0.1 * r2)
r = 0.1 * r2;
endif
endwhile
t = log (r);
y = lambda(k) .* r + log (sqrt (2*r.*((1-r) + r.*t)) ./ abs (r-1)) ./ t;
inv(k) = floor (y - 0.0218 ./ (y + 0.065 * lambda(k)));
endif
endfunction
%!shared x
%! x = [-1 0 0.5 1 2];
%!assert (poissinv (x, ones (1,5)), [NaN 0 1 Inf NaN])
%!assert (poissinv (x, 1), [NaN 0 1 Inf NaN])
%!assert (poissinv (x, [1 0 NaN 1 1]), [NaN NaN NaN Inf NaN])
%!assert (poissinv ([x(1:2) NaN x(4:5)], 1), [NaN 0 NaN Inf NaN])
## Test class of input preserved
%!assert (poissinv ([x, NaN], 1), [NaN 0 1 Inf NaN NaN])
%!assert (poissinv (single ([x, NaN]), 1), single ([NaN 0 1 Inf NaN NaN]))
%!assert (poissinv ([x, NaN], single (1)), single ([NaN 0 1 Inf NaN NaN]))
## Test input validation
%!error poissinv ()
%!error poissinv (1)
%!error poissinv (1,2,3)
%!error poissinv (ones (3), ones (2))
%!error poissinv (ones (2), ones (3))
%!error poissinv (i, 2)
%!error poissinv (2, i)