Mercurial > octave
annotate scripts/statistics/distributions/binoinv.m @ 23219:3ac9f9ecfae5 stable
maint: Update copyright dates.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 22 Feb 2017 12:39:29 -0500 |
parents | e9a0469dedd9 |
children | 092078913d54 |
rev | line source |
---|---|
23219
3ac9f9ecfae5
maint: Update copyright dates.
John W. Eaton <jwe@octave.org>
parents:
23083
diff
changeset
|
1 ## Copyright (C) 2016-2017 Lachlan Andrew |
22323
bac0d6f07a3e
maint: Update copyright notices for 2016.
John W. Eaton <jwe@octave.org>
parents:
21758
diff
changeset
|
2 ## Copyright (C) 2012-2016 Rik Wehbring |
20607
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
3 ## Copyright (C) 1995-2012 Kurt Hornik |
5410 | 4 ## |
5 ## This file is part of Octave. | |
6 ## | |
7 ## Octave is free software; you can redistribute it and/or modify it | |
8 ## under the terms of the GNU General Public License as published by | |
7016 | 9 ## the Free Software Foundation; either version 3 of the License, or (at |
10 ## your option) any later version. | |
5410 | 11 ## |
12 ## Octave is distributed in the hope that it will be useful, but | |
13 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
15 ## General Public License for more details. | |
16 ## | |
17 ## You should have received a copy of the GNU General Public License | |
7016 | 18 ## along with Octave; see the file COPYING. If not, see |
19 ## <http://www.gnu.org/licenses/>. | |
5410 | 20 |
21 ## -*- texinfo -*- | |
20852
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20609
diff
changeset
|
22 ## @deftypefn {} {} binoinv (@var{x}, @var{n}, @var{p}) |
20174
d9341b422488
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
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:
19833
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 | 27 ## @end deftypefn |
28 | |
5411 | 29 function inv = binoinv (x, n, p) |
5410 | 30 |
31 if (nargin != 3) | |
6046 | 32 print_usage (); |
5410 | 33 endif |
34 | |
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
35 if (! isscalar (n) || ! isscalar (p)) |
5410 | 36 [retval, x, n, p] = common_size (x, n, p); |
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 | 39 endif |
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 | 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 | 50 endif |
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))); |
20609
4e307c55a2b5
Use isempty () rather than any () for faster code in inverse statistical distributions.
Rik <rik@octave.org>
parents:
20607
diff
changeset
|
58 if (! isempty (k)) |
20607
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
59 x = x(k); |
5410 | 60 if (isscalar (n) && isscalar (p)) |
20607
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
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:
20174
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:
20174
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:
20174
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:
20174
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 |
20607
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
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:
20174
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:
20174
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:
20174
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:
20174
diff
changeset
|
71 endif |
5410 | 72 endif |
73 endif | |
74 | |
75 endfunction | |
13171
19b9f17d22af
Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
76 |
20607
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
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:
20174
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:
20174
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:
20174
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:
20174
diff
changeset
|
81 function [m, k] = scalar_binoinv (x, n, p) |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21031
diff
changeset
|
82 |
20607
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
83 k = 1:length (x); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
84 m = zeros (size (x)); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
85 prev_limit = 0; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
86 limit = 10; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
87 cdf = 0; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
88 v = 0; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
89 do |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
90 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:
20174
diff
changeset
|
91 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:
20174
diff
changeset
|
92 [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:
20174
diff
changeset
|
93 m(k) += prev_limit - 1; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
94 k = k(v == 0); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
95 |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
96 prev_limit = limit; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
97 limit += limit; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
98 until (isempty (k) || limit >= 1000) |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
99 |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
100 endfunction |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
101 |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
102 ## 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:
20174
diff
changeset
|
103 ## 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:
20174
diff
changeset
|
104 ## 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:
20174
diff
changeset
|
105 ## 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:
20174
diff
changeset
|
106 ## 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:
20174
diff
changeset
|
107 function [m, k] = vector_binoinv (x, n, p) |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21031
diff
changeset
|
108 |
20607
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
109 k = 1:length(x); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
110 m = zeros (size (x)); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
111 prev_limit = 0; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
112 limit = 10; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
113 cdf = 0; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
114 v = 0; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
115 do |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
116 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:
20174
diff
changeset
|
117 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:
20174
diff
changeset
|
118 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:
20174
diff
changeset
|
119 pdf = binopdf (xx, nn, pp); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
120 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:
20174
diff
changeset
|
121 cdf = cumsum (pdf, 2); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
122 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:
20174
diff
changeset
|
123 [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:
20174
diff
changeset
|
124 m(k) += prev_limit - 1; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
125 k = k(v == 0); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
126 |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
127 prev_limit = limit; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
128 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:
20174
diff
changeset
|
129 until (isempty (k) || limit >= 1000) |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
130 |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
131 endfunction |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
132 |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
133 ## Vectorized binary search. |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
134 ## 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:
20174
diff
changeset
|
135 ## answer is large. |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
136 ## 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:
20174
diff
changeset
|
137 ## 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:
20174
diff
changeset
|
138 ## 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:
20174
diff
changeset
|
139 function m = bin_search_binoinv (x, n, p) |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21031
diff
changeset
|
140 |
20607
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
141 k = 1:length (x); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
142 lower = zeros (size (x)); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
143 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:
20174
diff
changeset
|
144 while (any (k) && limit < 1e100) |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
145 cdf = binocdf (limit, n, p); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
146 k = (x > cdf); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
147 lower(k) = limit; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
148 limit += limit; |
21031
66a08c3cafe3
maint: Follow Octave coding conventions in m-files.
Rik <rik@octave.org>
parents:
20955
diff
changeset
|
149 endwhile |
20607
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
150 upper = max (2*lower, 1); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
151 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:
20174
diff
changeset
|
152 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:
20174
diff
changeset
|
153 mid = (upper + lower)/2; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
154 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:
20174
diff
changeset
|
155 r = (x <= cdf); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
156 upper(r) = mid(r); |
20955
77f5591878bf
maint: Use '! expr' rather than '!expr' to conform to coding guidelines.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
157 lower(! r) = mid(! r); |
20607
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
158 endfor |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
159 m = ceil (lower); |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
160 m(x > binocdf (m(:), n, p)) += 1; # fix off-by-one errors from binary search |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21031
diff
changeset
|
161 |
20607
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
162 endfunction |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
163 |
13171
19b9f17d22af
Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
164 |
19b9f17d22af
Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
165 %!shared x |
19b9f17d22af
Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
166 %! 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
|
167 %!assert (binoinv (x, 2*ones (1,5), 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
|
168 %!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
|
169 %!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
|
170 %!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
|
171 %!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
|
172 %!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
|
173 |
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
174 ## 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
|
175 %!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
|
176 %!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
|
177 %!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
|
178 %!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
|
179 |
20607
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
180 ## 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:
20174
diff
changeset
|
181 %!shared y, tol |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
182 %! y = magic (3) + 1; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
183 %! tol = 1; |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
184 %!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:
20174
diff
changeset
|
185 %!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:
20174
diff
changeset
|
186 %!assert (binoinv (binocdf (y, 2*y, 1./y), 2*y, 1./y), y, tol) |
9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
Lachlan Andrew <lachlanbis@gmail.com>
parents:
20174
diff
changeset
|
187 |
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
188 ## Test input validation |
13171
19b9f17d22af
Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
189 %!error binoinv () |
19b9f17d22af
Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
190 %!error binoinv (1) |
19b9f17d22af
Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
191 %!error binoinv (1,2) |
19b9f17d22af
Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
192 %!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
|
193 %!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
|
194 %!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
|
195 %!error binoinv (ones (2), ones (2), ones (3)) |
13171
19b9f17d22af
Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
196 %!error binoinv (i, 2, 2) |
19b9f17d22af
Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
197 %!error binoinv (2, i, 2) |
19b9f17d22af
Overhaul of statistical distribution functions
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
198 %!error binoinv (2, 2, i) |