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

annotate scripts/statistics/distributions/binoinv.m @ 20642:9d2023d1a63c

binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). * binoinv.m: Call new functions scalar_binoinv or vector_binoinv to calculate binoinv. If there are still uncalculated values then call bin_search_binoinv to perform binary search for remaining values. Add more BIST tests. * binoinv.m (scalar_binoinv): New subfunction to calculate binoinv for scalar x. Stops when x > 1000. * binoinv.m (vector_binoinv): New subfunction to calculate binoinv for scalar x. Stops when x > 1000.

author	Lachlan Andrew <lachlanbis@gmail.com>
date	Sun, 11 Oct 2015 19:49:40 -0700
parents	d9341b422488
children	4e307c55a2b5

rev	line source
20642 9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	1 ## Copyright (C) 2015 Lachlan Andrew
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	2 ## Copyright (C) 2012-2015 Rik Wehbring
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	3 ## Copyright (C) 1995-2012 Kurt Hornik
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	4 ##
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	5 ## This file is part of Octave.
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	6 ##
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	7 ## 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	8 ## 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	9 ## 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	10 ## your option) any later version.
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	11 ##
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	12 ## 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	13 ## WITHOUT ANY WARRANTY; without even the implied warranty of
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	14 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	15 ## General Public License for more details.
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	16 ##
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	17 ## 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	18 ## 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	19 ## <http://www.gnu.org/licenses/>.
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	20
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	21 ## -- texinfo --
5411 bee21f388110 [project @ 2005-07-13 17:53:44 by jwe] jwe parents: 5410 diff changeset	22 ## @deftypefn {Function File} {} binoinv (@var{x}, @var{n}, @var{p})
20209 d9341b422488 doc: Update more docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 19867 diff changeset	23 ## 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	24 ## at @var{x} of the binomial distribution with parameters
13949 0c15fece33ad doc: Clarify binomial and gamma distribution paramters Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13171 diff changeset	25 ## @var{n} and @var{p}, where @var{n} is the number of trials and
0c15fece33ad doc: Clarify binomial and gamma distribution paramters Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13171 diff changeset	26 ## @var{p} is the probability of success.
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	27 ## @end deftypefn
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	28
5411 bee21f388110 [project @ 2005-07-13 17:53:44 by jwe] jwe parents: 5410 diff changeset	29 function inv = binoinv (x, n, p)
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	30
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	31 if (nargin != 3)
6046 34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe] jwe parents: 5482 diff changeset	32 print_usage ();
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	33 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	34
19867 9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions. Rik <rik@octave.org> parents: 19731 diff changeset	35 if (! isscalar (n) \|\| ! isscalar (p))
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	36 [retval, x, n, p] = common_size (x, n, p);
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	37 if (retval > 0)
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	38 error ("binoinv: X, N, and P must be of common size or scalars");
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	39 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	40 endif
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	41
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	42 if (iscomplex (x) \|\| iscomplex (n) \|\| iscomplex (p))
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	43 error ("binoinv: X, N, and P must not be complex");
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	44 endif
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	45
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	46 if (isa (x, "single") \|\| isa (n, "single") \|\| isa (p, "single"));
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	47 inv = zeros (size (x), "single");
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	48 else
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	49 inv = zeros (size (x));
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	50 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	51
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	52 k = (!(x >= 0) \| !(x <= 1) \| !(n >= 0) \| (n != fix (n)) \|
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	53 !(p >= 0) \| !(p <= 1));
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	54 inv(k) = NaN;
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	55
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	56 k = find ((x >= 0) & (x <= 1) & (n >= 0) & (n == fix (n)
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	57 & (p >= 0) & (p <= 1)));
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	58 if (any (k))
20642 9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	59 x = x(k);
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	60 if (isscalar (n) && isscalar (p))
20642 9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	61 [inv(k), unfinished] = scalar_binoinv (x(:), n, p);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	62 k = k(unfinished);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	63 if (! isempty (k))
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	64 inv(k) = bin_search_binoinv (x(k), n, p);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	65 endif
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	66 else
20642 9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	67 [inv(k), unfinished] = vector_binoinv (x(:), n(:), p(:));
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	68 k = k(unfinished);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	69 if (! isempty (k))
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	70 inv(k) = bin_search_binoinv (x(k), n(k), p(k));
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	71 endif
5410 56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	72 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	73 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	74
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe] jwe parents: diff changeset	75 endfunction
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	76
20642 9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	77 ## Core algorithm to calculate the inverse binomial, for n and p real scalars
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	78 ## and y a column vector, and for which the output is not NaN or Inf.
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	79 ## Compute CDF in batches of doubling size until CDF > x, or answer > 500
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	80 ## Return the locations of unfinished cases in k.
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	81 function [m, k] = scalar_binoinv (x, n, p)
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	82 k = 1:length (x);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	83 m = zeros (size (x));
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	84 prev_limit = 0;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	85 limit = 10;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	86 cdf = 0;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	87 v = 0;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	88 do
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	89 cdf = binocdf (prev_limit:limit-1, n, p);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	90 r = bsxfun (@le, x(k), cdf);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	91 [v, m(k)] = max (r, [], 2); # find first instance of x <= cdf
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	92 m(k) += prev_limit - 1;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	93 k = k(v == 0);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	94
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	95 prev_limit = limit;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	96 limit += limit;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	97 until (isempty (k) \|\| limit >= 1000)
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	98
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	99 endfunction
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	100
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	101 ## Core algorithm to calculate the inverse binomial, for n, p, and y column
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	102 ## vectors, and for which the output is not NaN or Inf.
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	103 ## Compute CDF in batches of doubling size until CDF > x, or answer > 500
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	104 ## Return the locations of unfinished cases in k.
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	105 ## Calculates CDF by summing PDF, which is faster than calls to binocdf.
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	106 function [m, k] = vector_binoinv (x, n, p)
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	107 k = 1:length(x);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	108 m = zeros (size (x));
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	109 prev_limit = 0;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	110 limit = 10;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	111 cdf = 0;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	112 v = 0;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	113 do
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	114 xx = repmat (prev_limit:limit-1, [length(k), 1]);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	115 nn = kron (ones (1, limit-prev_limit), n(k));
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	116 pp = kron (ones (1, limit-prev_limit), p(k));
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	117 pdf = binopdf (xx, nn, pp);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	118 pdf(:,1) += cdf(v==0, end);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	119 cdf = cumsum (pdf, 2);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	120 r = bsxfun (@le, x(k), cdf);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	121 [v, m(k)] = max (r, [], 2); # find first instance of x <= cdf
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	122 m(k) += prev_limit - 1;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	123 k = k(v == 0);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	124
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	125 prev_limit = limit;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	126 limit += min (limit, max (1e4/numel (k), 10)); # limit memory use
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	127 until (isempty (k) \|\| limit >= 1000)
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	128
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	129 endfunction
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	130
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	131 ## Vectorized binary search.
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	132 ## Can handle vectors n and p, and is faster than the scalar case when the
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	133 ## answer is large.
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	134 ## Could be optimized to call binocdf only for a subset of the x at each stage,
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	135 ## but care must be taken to handle both scalar and vector n, p. Bookkeeping
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	136 ## may cost more than the extra computations.
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	137 function m = bin_search_binoinv (x, n, p)
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	138 k = 1:length (x);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	139 lower = zeros (size (x));
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	140 limit = 500; # lower bound on point at which prev phase finished
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	141 while (any (k) && limit < 1e100)
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	142 cdf = binocdf (limit, n, p);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	143 k = (x > cdf);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	144 lower(k) = limit;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	145 limit += limit;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	146 end
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	147 upper = max (2*lower, 1);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	148 k = find (lower != limit/2); # elements for which above loop finished
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	149 for i = 1:ceil (log2 (max (lower)))
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	150 mid = (upper + lower)/2;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	151 cdf = binocdf (floor(mid(:)), n, p);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	152 r = (x <= cdf);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	153 upper(r) = mid(r);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	154 lower(!r) = mid(!r);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	155 endfor
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	156 m = ceil (lower);
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	157 m(x > binocdf (m(:), n, p)) += 1; # fix off-by-one errors from binary search
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	158 endfunction
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	159
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	160
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	161 %!shared x
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	162 %! 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	163 %!assert (binoinv (x, 2ones (1,5), 0.5ones (1,5)), [NaN 0 1 2 NaN])
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	164 %!assert (binoinv (x, 2, 0.5*ones (1,5)), [NaN 0 1 2 NaN])
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	165 %!assert (binoinv (x, 2*ones (1,5), 0.5), [NaN 0 1 2 NaN])
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	166 %!assert (binoinv (x, 2*[0 -1 NaN 1.1 1], 0.5), [NaN NaN NaN NaN NaN])
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	167 %!assert (binoinv (x, 2, 0.5*[0 -1 NaN 3 1]), [NaN NaN NaN NaN NaN])
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	168 %!assert (binoinv ([x(1:2) NaN x(4:5)], 2, 0.5), [NaN 0 NaN 2 NaN])
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	169
19867 9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions. Rik <rik@octave.org> parents: 19731 diff changeset	170 ## 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	171 %!assert (binoinv ([x, NaN], 2, 0.5), [NaN 0 1 2 NaN NaN])
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	172 %!assert (binoinv (single ([x, NaN]), 2, 0.5), single ([NaN 0 1 2 NaN NaN]))
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	173 %!assert (binoinv ([x, NaN], single (2), 0.5), single ([NaN 0 1 2 NaN NaN]))
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	174 %!assert (binoinv ([x, NaN], 2, single (0.5)), single ([NaN 0 1 2 NaN NaN]))
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	175
20642 9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	176 ## Test accuracy, to within +/- 1 since it is a discrete distribution
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	177 %!shared y, tol
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	178 %! y = magic (3) + 1;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	179 %! tol = 1;
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	180 %!assert (binoinv (binocdf (1:10, 11, 0.1), 11, 0.1), 1:10, tol)
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	181 %!assert (binoinv (binocdf (1:10, 2(1:10), 0.1), 2(1:10), 0.1), 1:10, tol)
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	182 %!assert (binoinv (binocdf (y, 2y, 1./y), 2y, 1./y), y, tol)
9d2023d1a63c binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363). Lachlan Andrew <lachlanbis@gmail.com> parents: 20209 diff changeset	183
19867 9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions. Rik <rik@octave.org> parents: 19731 diff changeset	184 ## Test input validation
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	185 %!error binoinv ()
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	186 %!error binoinv (1)
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	187 %!error binoinv (1,2)
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	188 %!error binoinv (1,2,3,4)
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	189 %!error binoinv (ones (3), ones (2), ones (2))
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	190 %!error binoinv (ones (2), ones (3), ones (2))
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	191 %!error binoinv (ones (2), ones (2), ones (3))
13171 19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	192 %!error binoinv (i, 2, 2)
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	193 %!error binoinv (2, i, 2)
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	194 %!error binoinv (2, 2, i)
19b9f17d22af Overhaul of statistical distribution functions Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	195

Mercurial > octave-nkf

annotate scripts/statistics/distributions/binoinv.m @ 20642:9d2023d1a63c