time | Calls | line |
---|
| | 1 | function [x,xlo,xup] = gaminv(p,a,b,pcov,alpha)
|
| | 2 | %GAMINV Inverse of the gamma cumulative distribution function (cdf).
|
| | 3 | % X = GAMINV(P,A,B) returns the inverse cdf for the gamma distribution
|
| | 4 | % with shape A and scale B, evaluated at the values in P. The size of X
|
| | 5 | % is the common size of the input arguments. A scalar input functions as
|
| | 6 | % a constant matrix of the same size as the other inputs.
|
| | 7 | %
|
| | 8 | % [X,XLO,XUP] = GAMINV(P,A,B,PCOV,ALPHA) produces confidence bounds for
|
| | 9 | % X when the input parameters A and B are estimates. PCOV is a 2-by-2
|
| | 10 | % matrix containing the covariance matrix of the estimated parameters.
|
| | 11 | % ALPHA has a default value of 0.05, and specifies 100*(1-ALPHA)%
|
| | 12 | % confidence bounds. XLO and XUP are arrays of the same size as X
|
| | 13 | % containing the lower and upper confidence bounds.
|
| | 14 | %
|
| | 15 | % See also GAMCDF, GAMFIT, GAMLIKE, GAMPDF, GAMRND, GAMSTAT.
|
| | 16 |
|
| | 17 | % GAMINV uses Newton's method to find roots of GAMCDF(X,A,B) = P.
|
| | 18 |
|
| | 19 | % References:
|
| | 20 | % [1] Abramowitz, M. and Stegun, I.A. (1964) Handbook of Mathematical
|
| | 21 | % Functions, Dover, New York, section 26.1.
|
| | 22 | % [2] Evans, M., Hastings, N., and Peacock, B. (1993) Statistical
|
| | 23 | % Distributions, 2nd ed., Wiley.
|
| | 24 |
|
| | 25 | % Copyright 1993-2011 The MathWorks, Inc.
|
| | 26 |
|
| | 27 |
|
< 0.01 | 9 | 28 | if nargin < 2
|
| | 29 | error(message('stats:gaminv:TooFewInputs'));
|
< 0.01 | 9 | 30 | elseif nargin < 3
|
| | 31 | b = 1;
|
| | 32 | end
|
| | 33 |
|
| | 34 | % More checking if we need to compute confidence bounds.
|
< 0.01 | 9 | 35 | if nargout > 2
|
| | 36 | if nargin < 4
|
| | 37 | error(message('stats:gaminv:NeedCovariance'));
|
| | 38 | end
|
| | 39 | if ~isequal(size(pcov),[2 2])
|
| | 40 | error(message('stats:gaminv:BadCovarianceSize'));
|
| | 41 | end
|
| | 42 | if nargin < 5
|
| | 43 | alpha = 0.05;
|
| | 44 | elseif ~isnumeric(alpha) || numel(alpha) ~= 1 || alpha <= 0 || alpha >= 1
|
| | 45 | error(message('stats:gaminv:BadAlpha'));
|
| | 46 | end
|
| | 47 | end
|
| | 48 |
|
| | 49 | % Weed out any out of range parameters or edge/bad probabilities.
|
< 0.01 | 9 | 50 | try
|
< 0.01 | 9 | 51 | okAB = (0 < a & a < Inf) & (0 < b);
|
< 0.01 | 9 | 52 | k = (okAB & (0 < p & p < 1)); % GAMMAINCINV would handle 0 or 1, but the CI code won't
|
| | 53 | catch
|
| | 54 | error(message('stats:gaminv:InputSizeMismatch'));
|
| | 55 | end
|
< 0.01 | 9 | 56 | allOK = all(k(:));
|
| | 57 |
|
| | 58 | % Fill in NaNs for out of range cases, fill in edges cases when P is 0 or 1.
|
< 0.01 | 9 | 59 | if ~allOK
|
| | 60 | x = NaN(size(k),'like',internal.stats.dominantType(p,a,b)); % single if p, a, or b is
|
| | 61 | x(okAB & p == 0) = 0;
|
| | 62 | x(okAB & p == 1) = Inf;
|
| | 63 |
|
| | 64 | x(a==0) = 0;
|
| | 65 |
|
| | 66 | if nargout > 1
|
| | 67 | xlo = x; % NaNs or zeros or Infs
|
| | 68 | xup = x; % NaNs or zeros or Infs
|
| | 69 | end
|
| | 70 |
|
| | 71 | % Remove the bad/edge cases, leaving the easy cases. If there's
|
| | 72 | % nothing remaining, return.
|
| | 73 | if any(k(:))
|
| | 74 | if numel(p) > 1, p = p(k); end
|
| | 75 | if numel(a) > 1, a = a(k); end
|
| | 76 | if numel(b) > 1, b = b(k); end
|
| | 77 | else
|
| | 78 | return;
|
| | 79 | end
|
| | 80 | end
|
| | 81 |
|
| | 82 | % Call BETAINCINV to find a root of BETAINC(Q,A,B) = P
|
< 0.01 | 9 | 83 | q = gammaincinv(p,a);
|
| | 84 |
|
< 0.01 | 9 | 85 | tolerance = sqrt(eps(ones('like',q)));
|
< 0.01 | 9 | 86 | badcdf = ((abs(gammainc(q,a) - p)./p) > tolerance);
|
< 0.01 | 9 | 87 | if any(badcdf(:)) % cdf is too far off
|
| | 88 | didnt = find(badcdf, 1, 'first');
|
| | 89 | if numel(a) == 1, abad = a; else abad = a(didnt); end
|
| | 90 | if numel(b) == 1, bbad = b; else bbad = b(didnt); end
|
| | 91 | if numel(p) == 1, pbad = p; else pbad = p(didnt); end
|
| | 92 | warning(message('stats:gaminv:NoConvergence', sprintf( '%g', abad ), sprintf( '%g', bbad ), sprintf( '%g', pbad )));
|
| | 93 | end
|
| | 94 |
|
| | 95 | % Add in the scale factor, and broadcast the values to the correct place if
|
| | 96 | % need be.
|
< 0.01 | 9 | 97 | if allOK
|
< 0.01 | 9 | 98 | x = q .* b;
|
| | 99 | else
|
| | 100 | x(k) = q .* b;
|
| | 101 | end
|
| | 102 |
|
| | 103 | % Compute confidence bounds if requested.
|
< 0.01 | 9 | 104 | if nargout >= 2
|
| | 105 | logq = log(q);
|
| | 106 | dqda = -dgammainc(q,a) ./ exp((a-1).*logq - q - gammaln(a));
|
| | 107 |
|
| | 108 | % Approximate the variance of x=q*b on the log scale.
|
| | 109 | % dlogx/da = dlogx/dq * dqda = dqda/q
|
| | 110 | % dlogx/db = 1/b
|
| | 111 | logx = logq + log(b);
|
| | 112 | varlogx = pcov(1,1).*(dqda./q).^2 + 2.*pcov(1,2).*dqda./(b.*q) + pcov(2,2)./(b.^2);
|
| | 113 | if any(varlogx(:) < 0)
|
| | 114 | error(message('stats:gaminv:BadCovariancePosDef'));
|
| | 115 | end
|
| | 116 | z = -norminv(alpha/2);
|
| | 117 | halfwidth = z * sqrt(varlogx);
|
| | 118 |
|
| | 119 | % Convert back to original scale
|
| | 120 | if allOK
|
| | 121 | xlo = exp(logx - halfwidth);
|
| | 122 | xup = exp(logx + halfwidth);
|
| | 123 | else
|
| | 124 | xlo(k) = exp(logx - halfwidth);
|
| | 125 | xup(k) = exp(logx + halfwidth);
|
| | 126 | end
|
| | 127 | end
|