annotate scripts/statistics/distributions/nbininv.m @ 20643:d6d04088ac9e

nbininv.m: Increase speed (85X) and accuracy of function (bug #34363). * nbininv.m: Call new function scalar_nbininv to calculate nbininv for scalar. If there are still uncalculated values then call bin_search_nbininv. Call bin_search_nbininv directly for vectors. Add more BIST tests. * nbininv.m (scalar_binoinv): New subfunction to calculate nbininv for scalar x. Stops when x > 1000. * nbininv.m (bin_search_nbininv): New subfunction to do binary search for nbininv.
author Lachlan Andrew <lachlanbis@gmail.com>
date Sun, 11 Oct 2015 20:33:37 -0700
parents d9341b422488
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
20643
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
1 ## Copyright (C) 2015 Lachlan Andrew
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
2 ## Copyright (C) 2012-2015 Rik Wehbring
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
3 ## Copyright (C) 1995-2012 Kurt Hornik
6349
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
4 ##
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
5 ## This file is part of Octave.
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
6 ##
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
7 ## Octave is free software; you can redistribute it and/or modify it
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
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: 6349
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: 6349
diff changeset
10 ## your option) any later version.
6349
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
11 ##
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
12 ## Octave is distributed in the hope that it will be useful, but
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
13 ## WITHOUT ANY WARRANTY; without even the implied warranty of
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
14 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
15 ## General Public License for more details.
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
16 ##
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
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: 6349
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: 6349
diff changeset
19 ## <http://www.gnu.org/licenses/>.
6349
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
20
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
21 ## -*- texinfo -*-
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
22 ## @deftypefn {Function File} {} nbininv (@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 negative binomial distribution with parameters
d9341b422488 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 19867
diff changeset
25 ## @var{n} and @var{p}.
6349
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
26 ##
20209
d9341b422488 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 19867
diff changeset
27 ## When @var{n} is integer this is the Pascal distribution.
d9341b422488 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 19867
diff changeset
28 ## When @var{n} is extended to real numbers this is the Polya distribution.
19627
446c46af4b42 strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents: 17744
diff changeset
29 ##
20209
d9341b422488 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 19867
diff changeset
30 ## The number of failures in a Bernoulli experiment with success probability
d9341b422488 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 19867
diff changeset
31 ## @var{p} before the @var{n}-th success follows this distribution.
6349
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
32 ## @end deftypefn
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
33
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
34 function inv = nbininv (x, n, p)
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
35
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
36 if (nargin != 3)
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
37 print_usage ();
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
38 endif
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
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 (n) || ! isscalar (p))
6349
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
41 [retval, x, n, p] = common_size (x, n, p);
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
42 if (retval > 0)
13171
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
43 error ("nbininv: X, N, and P must be of common size or scalars");
6349
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
44 endif
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
45 endif
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
46
13171
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
47 if (iscomplex (x) || iscomplex (n) || iscomplex (p))
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
48 error ("nbininv: X, N, and P must not be complex");
6349
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
49 endif
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
50
13171
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
51 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
52 inv = zeros (size (x), "single");
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
53 else
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
54 inv = zeros (size (x));
6349
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
55 endif
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
56
13171
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
57 k = (isnan (x) | (x < 0) | (x > 1) | isnan (n) | (n < 1) | (n == Inf)
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
58 | isnan (p) | (p < 0) | (p > 1));
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
59 inv(k) = NaN;
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
60
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
61 k = (x == 1) & (n > 0) & (n < Inf) & (p >= 0) & (p <= 1);
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
62 inv(k) = Inf;
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
63
6349
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
64 k = find ((x >= 0) & (x < 1) & (n > 0) & (n < Inf)
13171
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
65 & (p > 0) & (p <= 1));
20643
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
66 if (! isempty (k))
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
67 x = x(k);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
68 m = zeros (size (k));
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
69 if (isscalar (n) && isscalar (p))
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
70 [m, unfinished] = scalar_nbininv (x(:), n, p);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
71 m(unfinished) = bin_search_nbininv (x(unfinished), n, p);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
72 else
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
73 m = bin_search_nbininv (x, n(k), p(k));
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
74 endif
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
75 inv(k) = m;
6349
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
76 endif
20643
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
77
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
78 endfunction
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
79
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
80
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
81 ## Core algorithm to calculate the inverse negative binomial, for n and p real
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
82 ## scalars and y a column vector, and for which the output is not NaN or Inf.
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
83 ## Compute CDF in batches of doubling size until CDF > x, or answer > 500.
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
84 ## Return the locations of unfinished cases in k.
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
85 function [m, k] = scalar_nbininv (x, n, p)
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
86 k = 1:length (x);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
87 m = zeros (size (x));
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
88 prev_limit = 0;
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
89 limit = 10;
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
90 do
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
91 cdf = nbincdf (prev_limit:limit, n, p);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
92 r = bsxfun (@le, x(k), cdf);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
93 [v, m(k)] = max (r, [], 2); # find first instance of x <= cdf
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
94 m(k) += prev_limit - 1;
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
95 k = k(v == 0);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
96
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
97 prev_limit = limit;
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
98 limit += limit;
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
99 until (isempty (k) || limit >= 1000)
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
100
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
101 endfunction
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
102
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
103 ## Vectorized binary search.
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
104 ## Can handle vectors n and p, and is faster than the scalar case when the
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
105 ## answer is large.
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
106 ## Could be optimized to call nbincdf only for a subset of the x at each stage,
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
107 ## but care must be taken to handle both scalar and vector n,p. Bookkeeping
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
108 ## may cost more than the extra computations.
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
109 function m = bin_search_nbininv (x, n, p)
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
110 k = 1:length (x);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
111 lower = zeros (size (x));
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
112 limit = 1;
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
113 while (any (k) && limit < 1e100)
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
114 cdf = nbincdf (limit, n, p);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
115 k = (x > cdf);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
116 lower(k) = limit;
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
117 limit += limit;
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
118 end
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
119 upper = max (2*lower, 1);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
120 k = find (lower != limit/2); # elements for which above loop finished
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
121 for i = 1:ceil (log2 (max (lower)))
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
122 mid = (upper + lower)/2;
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
123 cdf = nbincdf (floor (mid), n, p);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
124 r = (x <= cdf);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
125 upper(r) = mid(r);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
126 lower(!r) = mid(!r);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
127 endfor
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
128 m = ceil (lower);
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
129 m(x > nbincdf (m, n, p)) += 1; # fix off-by-one errors from binary search
6349
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
130
859f7aaea254 [project @ 2007-02-23 14:20:42 by dbateman]
dbateman
parents:
diff changeset
131 endfunction
13171
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
132
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
133
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
134 %!shared x
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
135 %! x = [-1 0 3/4 1 2];
14363
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14138
diff changeset
136 %!assert (nbininv (x, ones (1,5), 0.5*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
137 %!assert (nbininv (x, 1, 0.5*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
138 %!assert (nbininv (x, ones (1,5), 0.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
139 %!assert (nbininv (x, [1 0 NaN Inf 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
140 %!assert (nbininv (x, [1 0 1.5 Inf 1], 0.5), [NaN NaN 2 NaN NaN])
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14138
diff changeset
141 %!assert (nbininv (x, 1, 0.5*[1 -Inf NaN Inf 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
142 %!assert (nbininv ([x(1:2) NaN x(4:5)], 1, 0.5), [NaN 0 NaN Inf NaN])
13171
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
143
19867
9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19731
diff changeset
144 ## 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
145 %!assert (nbininv ([x, NaN], 1, 0.5), [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
146 %!assert (nbininv (single ([x, NaN]), 1, 0.5), 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
147 %!assert (nbininv ([x, NaN], single (1), 0.5), 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
148 %!assert (nbininv ([x, NaN], 1, single (0.5)), single ([NaN 0 1 Inf NaN NaN]))
13171
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
149
20643
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
150 ## Test accuracy, to within +/- 1 since it is a discrete distribution
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
151 %!shared y, tol
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
152 %! y = magic (3) + 1;
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
153 %! tol = 1;
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
154 %!assert (nbininv (nbincdf (1:10, 3, 0.1), 3, 0.1), 1:10, tol)
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
155 %!assert (nbininv (nbincdf (1:10, 3./(1:10), 0.1), 3./(1:10), 0.1), 1:10, tol)
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
156 %!assert (nbininv (nbincdf (y, 3./y, 1./y), 3./y, 1./y), y, tol)
d6d04088ac9e nbininv.m: Increase speed (85X) and accuracy of function (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents: 20209
diff changeset
157
19867
9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19731
diff changeset
158 ## Test input validation
13171
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
159 %!error nbininv ()
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
160 %!error nbininv (1)
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
161 %!error nbininv (1,2)
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
162 %!error nbininv (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
163 %!error nbininv (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
164 %!error nbininv (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
165 %!error nbininv (ones (2), ones (2), ones (3))
13171
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
166 %!error nbininv (i, 2, 2)
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
167 %!error nbininv (2, i, 2)
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
168 %!error nbininv (2, 2, i)
19b9f17d22af Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
169