This is a static copy of a profile report

Home

gaminv (Calls: 9, Time: 0.007 sec)
Generated 28-May-2016 14:58:47 using performance time.
function in file /home/johs/MATLAB/R2015b/toolbox/stats/stats/gaminv.m
Copy to new window for comparing multiple runs

Parents (calling functions)

Function NameFunction TypeCalls
chi2invfunction9
Lines where the most time was spent

Line NumberCodeCallsTotal Time% TimeTime Plot
86
badcdf = ((abs(gammainc(q,a) -...
90.002 s35.3%
83
q = gammaincinv(p,a);
90.002 s21.6%
50
try
90.001 s8.8%
52
k = (okAB & (0 < p &...
90.001 s8.4%
97
if allOK
90.000 s7.0%
All other lines  0.001 s18.9%
Totals  0.007 s100% 
Children (called functions)
No children
Code Analyzer results
No Code Analyzer messages.
Coverage results
Show coverage for parent directory
Total lines in function127
Non-code lines (comments, blank lines)52
Code lines (lines that can run)75
Code lines that did run15
Code lines that did not run60
Coverage (did run/can run)20.00 %
Function listing
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