octave-nkf: scripts/statistics/distributions/poissinv.m annotate

annotate scripts/statistics/distributions/poissinv.m @ 20644:4e307c55a2b5

Use isempty () rather than any () for faster code in inverse statistical distributions. betainv.m, binoinv.m, gaminv.m, poissinv.m: Use '! isempty (k)' rather than 'any (k)' for faster code.

author	Rik <rik@octave.org>
date	Sun, 11 Oct 2015 21:09:41 -0700
parents	b3f56ed8d1f9
children

rev	line source
20485 b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	1 ## Copyright (C) 1995-2015 Kurt Hornik
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	2 ## Copyright (C) 2015 Lachlan Andrew
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	3 ## Copyright (C) 2014 Mike Giles
14138 72c96de7a403 maint: update copyright notices for 2012 John W. Eaton <jwe@octave.org> parents: 13171 diff changeset	4 ## Copyright (C) 2012 Rik Wehbring
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	5 ##
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	6 ## This file is part of Octave.
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	7 ##
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	8 ## Octave is free software; you can redistribute it and/or modify it
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	9 ## under the terms of the GNU General Public License as published by
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6046 diff changeset	10 ## the Free Software Foundation; either version 3 of the License, or (at
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6046 diff changeset	11 ## your option) any later version.
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	12 ##
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	13 ## Octave is distributed in the hope that it will be useful, but
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	14 ## WITHOUT ANY WARRANTY; without even the implied warranty of
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	15 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	16 ## General Public License for more details.
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	17 ##
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	18 ## You should have received a copy of the GNU General Public License
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6046 diff changeset	19 ## along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6046 diff changeset	20 ## <http://www.gnu.org/licenses/>.
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	21
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	22 ## -- texinfo --
5411 bee21f388110 [project @ 2005-07-13 17:53:44 by jwe] jwe parents: 5410 diff changeset	23 ## @deftypefn {Function File} {} poissinv (@var{x}, @var{lambda})
20209 d9341b422488 doc: Update more docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 19867 diff changeset	24 ## For each element of @var{x}, compute the quantile (the inverse of the CDF)
d9341b422488 doc: Update more docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 19867 diff changeset	25 ## at @var{x} of the Poisson distribution with parameter @var{lambda}.
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	26 ## @end deftypefn
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	27
20485 b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	28 ## Author: Lachlan <lachlanbis@gmail.com>
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	29 ## based on code by
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	30 ## KH <Kurt.Hornik@wu-wien.ac.at>
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	31 ## Mike Giles <mike.giles@maths.ox.ac.uk>
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	32 ## Description: Quantile function of the Poisson distribution
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	33
11469 c776f063fefe Overhaul m-script files to use common variable name between code and documentation. Rik <octave@nomad.inbox5.com> parents: 10549 diff changeset	34 function inv = poissinv (x, lambda)
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	35
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	36 if (nargin != 2)
6046 34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe] jwe parents: 5482 diff changeset	37 print_usage ();
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	38 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	39
19867 9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions. Rik <rik@octave.org> parents: 19731 diff changeset	40 if (! isscalar (lambda))
11469 c776f063fefe Overhaul m-script files to use common variable name between code and documentation. Rik <octave@nomad.inbox5.com> parents: 10549 diff changeset	41 [retval, x, lambda] = common_size (x, lambda);
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	42 if (retval > 0)
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	43 error ("poissinv: X and LAMBDA must be of common size or scalars");
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	44 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	45 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	46
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	47 if (iscomplex (x) \|\| iscomplex (lambda))
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	48 error ("poissinv: X and LAMBDA must not be complex");
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	49 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	50
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	51 if (isa (x, "single") \|\| isa (lambda, "single"))
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	52 inv = zeros (size (x), "single");
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	53 else
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	54 inv = zeros (size (x));
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	55 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	56
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	57 k = (x < 0) \| (x > 1) \| isnan (x) \| !(lambda > 0);
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	58 inv(k) = NaN;
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	59
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	60 k = (x == 1) & (lambda > 0);
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	61 inv(k) = Inf;
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	62
20485 b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	63 k = (x > 0) & (x < 1) & (lambda > 0);
20644 4e307c55a2b5 Use isempty () rather than any () for faster code in inverse statistical distributions. Rik <rik@octave.org> parents: 20485 diff changeset	64 if (any (k(:)))
4e307c55a2b5 Use isempty () rather than any () for faster code in inverse statistical distributions. Rik <rik@octave.org> parents: 20485 diff changeset	65 limit = 20; # After 'limit' iterations, use approx
20485 b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	66 if (isscalar (lambda))
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	67 cdf = [(cumsum (poisspdf (0:limit-1,lambda))), 2];
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	68 y = x(:); # force to column
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	69 r = bsxfun (@le, y(k), cdf);
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	70 [~, inv(k)] = max (r, [], 2); # find first instance of x <= cdf
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	71 inv(k) -= 1;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	72 else
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	73 kk = find (k);
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	74 cdf = exp (-lambda(kk));
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	75 for i = 1:limit
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	76 m = find (cdf < x(kk));
20644 4e307c55a2b5 Use isempty () rather than any () for faster code in inverse statistical distributions. Rik <rik@octave.org> parents: 20485 diff changeset	77 if (isempty (m))
4e307c55a2b5 Use isempty () rather than any () for faster code in inverse statistical distributions. Rik <rik@octave.org> parents: 20485 diff changeset	78 break;
4e307c55a2b5 Use isempty () rather than any () for faster code in inverse statistical distributions. Rik <rik@octave.org> parents: 20485 diff changeset	79 else
20485 b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	80 inv(kk(m)) += 1;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	81 cdf(m) += poisspdf (i, lambda(kk(m)));
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	82 endif
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	83 endfor
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	84 endif
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	85
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	86 ## Use Mike Giles's magic when inv isn't < limit
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	87 k &= (inv == limit);
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	88 if (any (k(:)))
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	89 if (isscalar (lambda))
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	90 lam = repmat (lambda, size (x));
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	91 else
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	92 lam = lambda;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	93 endif
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	94 inv(k) = analytic_approx (x(k), lam(k));
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	95 endif
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	96 endif
19627 446c46af4b42 strip trailing whitespace from most source files John W. Eaton <jwe@octave.org> parents: 17744 diff changeset	97
20485 b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	98 endfunction
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	99
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	100
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	101 ## The following is based on Mike Giles's CUDA implementation,
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	102 ## [http://people.maths.ox.ac.uk/gilesm/codes/poissinv/poissinv_cuda.h]
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	103 ## which is copyright by the University of Oxford
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	104 ## and is provided under the terms of the GNU GPLv3 license:
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	105 ## http://www.gnu.org/licenses/gpl.html
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	106
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	107 function inv = analytic_approx (x, lambda)
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	108 s = norminv (x, 0, 1) ./ sqrt (lambda);
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	109 k = (s > -0.6833501) & (s < 1.777993);
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	110 ## use polynomial approximations in central region
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	111 if (any (k))
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	112 lam = lambda(k);
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	113 if (isscalar (s))
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	114 sk = s;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	115 else
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	116 sk = s(k);
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	117 endif
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	118
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	119 ## polynomial approximation to f^{-1}(s) - 1
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	120 rm = 2.82298751e-07;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	121 rm = -2.58136133e-06 + rm.*sk;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	122 rm = 1.02118025e-05 + rm.*sk;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	123 rm = -2.37996199e-05 + rm.*sk;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	124 rm = 4.05347462e-05 + rm.*sk;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	125 rm = -6.63730967e-05 + rm.*sk;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	126 rm = 0.000124762566 + rm.*sk;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	127 rm = -0.000256970731 + rm.*sk;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	128 rm = 0.000558953132 + rm.*sk;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	129 rm = -0.00133129194 + rm.*sk;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	130 rm = 0.00370367937 + rm.*sk;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	131 rm = -0.0138888706 + rm.*sk;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	132 rm = 0.166666667 + rm.*sk;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	133 rm = sk + sk.(rm.sk);
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	134
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	135 ## polynomial approximation to correction c0(r)
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	136
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	137 t = 1.86386867e-05;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	138 t = -0.000207319499 + t.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	139 t = 0.0009689451 + t.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	140 t = -0.00247340054 + t.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	141 t = 0.00379952985 + t.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	142 t = -0.00386717047 + t.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	143 t = 0.00346960934 + t.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	144 t = -0.00414125511 + t.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	145 t = 0.00586752093 + t.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	146 t = -0.00838583787 + t.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	147 t = 0.0132793933 + t.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	148 t = -0.027775536 + t.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	149 t = 0.333333333 + t.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	150
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	151 ## O(1/lam) correction
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	152
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	153 y = -0.00014585224;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	154 y = 0.00146121529 + y.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	155 y = -0.00610328845 + y.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	156 y = 0.0138117964 + y.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	157 y = -0.0186988746 + y.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	158 y = 0.0168155118 + y.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	159 y = -0.013394797 + y.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	160 y = 0.0135698573 + y.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	161 y = -0.0155377333 + y.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	162 y = 0.0174065334 + y.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	163 y = -0.0198011178 + y.*rm;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	164 y ./= lam;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	165
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	166 inv(k) = floor (lam + (y+t)+lam.*rm);
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	167 endif
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	168
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	169 k = !k & (s > -sqrt (2));
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	170 if (any (k))
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	171 ## Newton iteration
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	172 r = 1 + s(k);
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	173 r2 = r + 1;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	174 while (any (abs (r - r2) > 1e-5))
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	175 t = log (r);
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	176 r2 = r;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	177 s2 = sqrt (2 * ((1-r) + r.*t));
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	178 s2(r<1) *= -1;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	179 r = r2 - (s2 - s(k)) .* s2 ./ t;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	180 if (r < 0.1 * r2)
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	181 r = 0.1 * r2;
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	182 endif
20485 b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	183 endwhile
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	184 t = log (r);
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	185 y = lambda(k) .* r + log (sqrt (2r.((1-r) + r.*t)) ./ abs (r-1)) ./ t;
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	186 inv(k) = floor (y - 0.0218 ./ (y + 0.065 * lambda(k)));
b3f56ed8d1f9 poissinv.m: Overhaul function for accuracy and performance (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	187 endif
20644 4e307c55a2b5 Use isempty () rather than any () for faster code in inverse statistical distributions. Rik <rik@octave.org> parents: 20485 diff changeset	188
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	189 endfunction
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	190
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	191
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	192 %!shared x
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	193 %! x = [-1 0 0.5 1 2];
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	194 %!assert (poissinv (x, ones (1,5)), [NaN 0 1 Inf NaN])
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	195 %!assert (poissinv (x, 1), [NaN 0 1 Inf NaN])
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	196 %!assert (poissinv (x, [1 0 NaN 1 1]), [NaN NaN NaN Inf NaN])
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	197 %!assert (poissinv ([x(1:2) NaN x(4:5)], 1), [NaN 0 NaN Inf NaN])
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	198
19867 9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions. Rik <rik@octave.org> parents: 19731 diff changeset	199 ## Test class of input preserved
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	200 %!assert (poissinv ([x, NaN], 1), [NaN 0 1 Inf NaN NaN])
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	201 %!assert (poissinv (single ([x, NaN]), 1), single ([NaN 0 1 Inf NaN NaN]))
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	202 %!assert (poissinv ([x, NaN], single (1)), single ([NaN 0 1 Inf NaN NaN]))
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	203
19867 9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions. Rik <rik@octave.org> parents: 19731 diff changeset	204 ## Test input validation
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	205 %!error poissinv ()
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	206 %!error poissinv (1)
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	207 %!error poissinv (1,2,3)
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	208 %!error poissinv (ones (3), ones (2))
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	209 %!error poissinv (ones (2), ones (3))
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	210 %!error poissinv (i, 2)
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	211 %!error poissinv (2, i)
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	212

Mercurial > octave-nkf

annotate scripts/statistics/distributions/poissinv.m @ 20644:4e307c55a2b5