Mercurial > octave-nkf

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use isempty () rather than any () for faster code in inverse statistical distributions. betainv.m, binoinv.m, gaminv.m, poissinv.m: Use '! isempty (k)' rather than 'any (k)' for faster code.
author Rik <rik@octave.org>
date Sun, 11 Oct 2015 21:09:41 -0700
parents b3f56ed8d1f9
children
comparison
equal deleted inserted replaced
20643:d6d04088ac9e 20644:4e307c55a2b5
59 59
60 k = (x == 1) & (lambda > 0); 60 k = (x == 1) & (lambda > 0);
61 inv(k) = Inf; 61 inv(k) = Inf;
62 62
63 k = (x > 0) & (x < 1) & (lambda > 0); 63 k = (x > 0) & (x < 1) & (lambda > 0);
64 64 if (any (k(:)))
65 limit = 20; # After 'limit' iterations, use approx 65 limit = 20; # After 'limit' iterations, use approx
66 if (! isempty (k))
67 if (isscalar (lambda)) 66 if (isscalar (lambda))
68 cdf = [(cumsum (poisspdf (0:limit-1,lambda))), 2]; 67 cdf = [(cumsum (poisspdf (0:limit-1,lambda))), 2];
69 y = x(:); # force to column 68 y = x(:); # force to column
70 r = bsxfun (@le, y(k), cdf); 69 r = bsxfun (@le, y(k), cdf);
71 [~, inv(k)] = max (r, [], 2); # find first instance of x <= cdf 70 [~, inv(k)] = max (r, [], 2); # find first instance of x <= cdf
73 else 72 else
74 kk = find (k); 73 kk = find (k);
75 cdf = exp (-lambda(kk)); 74 cdf = exp (-lambda(kk));
76 for i = 1:limit 75 for i = 1:limit
77 m = find (cdf < x(kk)); 76 m = find (cdf < x(kk));
78 if (any (m)) 77 if (isempty (m))
78 break;
79 else
79 inv(kk(m)) += 1; 80 inv(kk(m)) += 1;
80 cdf(m) += poisspdf (i, lambda(kk(m))); 81 cdf(m) += poisspdf (i, lambda(kk(m)));
81 else
82 break;
83 endif 82 endif
84 endfor 83 endfor
85 endif 84 endif
86 85
87 ## Use Mike Giles's magic when inv isn't < limit 86 ## Use Mike Giles's magic when inv isn't < limit
184 endwhile 183 endwhile
185 t = log (r); 184 t = log (r);
186 y = lambda(k) .* r + log (sqrt (2*r.*((1-r) + r.*t)) ./ abs (r-1)) ./ t; 185 y = lambda(k) .* r + log (sqrt (2*r.*((1-r) + r.*t)) ./ abs (r-1)) ./ t;
187 inv(k) = floor (y - 0.0218 ./ (y + 0.065 * lambda(k))); 186 inv(k) = floor (y - 0.0218 ./ (y + 0.065 * lambda(k)));
188 endif 187 endif
188
189 endfunction 189 endfunction
190 190
191 191
192 %!shared x 192 %!shared x
193 %! x = [-1 0 0.5 1 2]; 193 %! x = [-1 0 0.5 1 2];