Compute the Inverse Gamma PDF, CDF, and iCDF in MATLAB Using Built-In Functions for the Gamma Distribution

I wrote about computing the inverse Gamma PDF and CDF in MATLAB using the known formula. For something I am working on, I have to compute the inverse CDF (iCDF) for the inverse Gamma distribution, which is not an easy task. As a result, I was thinking about using the built-in gampdf, gamcdf, and gaminv functions to compute the inverse Gamma PDF, CDF, and iCDF.

PDF

Given a random variable (RV) Y ~ Gamma(a, 1/b), define the transformation

$g(Y) = 1/Y$

and let

$X = g(Y) = 1/Y$.

Based on the definitions, X ~ Inv-Gamma(a, b). Because the transformation is invertible, the PDF for X, denoted by

$P_{X}(x)$ 

can be represented using the PDF for Y, denoted by 

$P_{Y}(y)$: 

$P_{X}(x) = \frac1{|dx/dy|}P_{Y}(y)\large|_{y = g^{-1}(x)}$. 

This leads to

$P_{X}(x) = y^2P_{Y}(y)\large|_{y = g^{-1}(x)} = P_{Y}(\frac1x)/x^2$. 

Therefore, in MATLAB, the inverse Gamma PDF for x for a shape parameter a and scale parameter b can be computed using gampdf(y,a,1/b)./(x.^2), or gampdf(1./x,a,1/b)/(x.^2).

A function can be created for this so that the similar code does not have to be rewritten every time when computing the PDF:

function [ Y ] = inversegampdfgam( X,A,B )
%inversegampdfgam Inverse gamma probability density function.
%   Y = inversegampdfgam(X,A,B) returns the inverse gamma probability
%   density function with shape and scale parameters A and B, respectively,
%   at the values in X.

%Y = B^A/gamma(A)*X.^(-A-1).*exp(-B./X);
Y = gampdf(1./X,A,1/B)./(X.^2);

end

Compare the results with those from computing the inverse Gamma PDF directly by finding the absolute difference:

In the figure above, the vertical axis represents the absolute difference between the results obtained using the formula and those obtained using gampdf. It can be seen that the difference is small (possibly resulting from estimation) and negligible.

CDF

Given a random variable (RV) Y ~ Gamma(a, 1/b), let

$X = 1/Y$.

Based on the definitions, X ~ Inv-Gamma(a, b). 

For X ≤ x, it follows that 1/Y ≤ 1/y, which leads to y ≤ Y. As a result, the probability for X ≤ x is equal to that of y ≤ Y. Therefore, the CDF for an inverse Gamma distribution can be computed using the iCDF for a Gamma distribution. In MATLAB, the inverse Gamma CDF for x for a shape parameter a and scale parameter b can then be computed using 1 - gamcdf(y,a,1/b), or 1 - gamcdf(1./x,a,1/b).

A function can be created for this so that the similar code does not have to be rewritten every time when computing the CDF:

function [ P ] = inversegamcdfgam( X,A,B )
%inversegamcdfgam Inverse gamma cumulative distribution function.
%   Y = inversegamcdfgam(X,A,B) returns the inverse gamma cumulative
%   distribution function with shape and scale parameters A and B,
%   respectively, at the values in X. The size of P is the common size of
%   the input arguments. A scalar input functions is a constant matrix of
%   the same size as the other inputs.

%P = gammainc(B./X,A,'upper');
P = 1 - gamcdf(1./X,A,1/B);

end

Compare the results with those from computing the inverse Gamma CDF directly by finding the absolute difference:

In the figure above, the vertical axis represents the absolute difference between the results obtained using the formula and those obtained using gamcdf. It can be seen that the difference is small (possibly resulting from estimation) and negligible.

iCDF

Based on the results for the CDF, the iCDF for the inverse Gamma distribution can be computed using the iCDF for the Gamma distribution.

Given a random variable (RV) Y ~ Gamma(a, 1/b), let

$X = g(Y) = 1/Y$.

Based on the definitions, X ~ Inv-Gamma(a, b).

Denote the CDF for X by F(x|a, b) and that for Y by F(y|a, 1/b), a being the shape parameter for X and b the scale parameter. It is known that

$F(x|a, b) = 1 - F(y|a, 1/b)$. 

Based on this, given 1 - F(x|a, b), the iCDF value for F(y|a, 1/b) is y, which is 1/x. Therefore, given the CDF value p, shape parameter a, and scale parameter b, the corresponding inverse Gamma RV x can be found in MATLAB using 1./gaminv(1-p,a,1/b).

A function can be created for this so that the similar code does not have to be rewritten every time when computing the PDF:

function [ X ] = inversegaminv( P,A,B )
%inversegaminv Inverse of the inverse gamma cumulative distribution
%function (cdf).
%   X = inversegaminv(P,A,B) returns the inverse cdf for the inverse gamma
%   distribution with shape A and scale B, evaluated at the values in P.
%   The size of X is the common size of the input arguments. A scalar input
%   functions is a constant matrix of the same size as the other inputs.

X = 1./gaminv(1 - P,A,1./B);

end

Given a set of CDF values, compare the CDF values for the iCDF results:

In the figure above, the vertical axis represents the absolute difference between the given CDF values and the CDF values computed using the formula for the iCDF results for the given CDF values. It can be seen that the difference is small (possibly resulting from estimation) and negligible.

-----

The above has not been peer-reviewed, but I will continue the project I am working on using the above information for now. If it is wrong, I will definitely see something wrong for the project I am working on (this is not a project that will cause any harm, so it is fine for me to continue for now).

Any constructive feedback will be appreciated.