GCIN
Evaluates the inverse of a general continuous cumulative distribution function given ordinates of the density.
Required Arguments
P — Probability for which the inverse of the distribution function is to be evaluated. (Input)
P must be in the open interval (0.0, 1.0).
X — Array containing the abscissas or the endpoints. (Input)
If IOPT = 1 or 3, X is of length 2. If IOPT = 2 or 4, X is of length M. For IOPT = 1 or 3, X(1) contains the lower endpoint of the support of the distribution and X(2) is the upper endpoint. For IOPT = 2 or 4, X contains, in strictly increasing order, the abscissas such that X(I) corresponds to F(I).
F — Vector of length M containing the probability density ordinates corresponding to increasing abscissas. (Input)
If IOPT = 1 or 3, for I = 1, 2, …, M, F(I) corresponds to X(1) + (I ‑ 1) * (X(2) ‑ X(1))/(M ‑ 1); otherwise, F and X correspond one for one.
GCIN — Function value. (Output)
The probability that a random variable whose density is given in F takes a value less than or equal to GCIN is P.
Optional Arguments
IOPT — Indicator of the method of interpolation. (Input)
Default: IOPT = 1.
IOPT | Interpolation Method |
|---|
1 | Linear interpolation with equally spaced abscissas. |
2 | Linear interpolation with possibly unequally spaced abscissas. |
3 | A cubic spline is fitted to equally spaced abscissas. |
4 | A cubic spline is fitted to possibly unequally spaced abscissas. |
M — Number of ordinates of the density supplied. (Input)
M must be greater than 1 for linear interpolation (IOPT = 1 or 2) and greater than 3 if a curve is fitted through the ordinates (IOPT = 3 or 4).
Default: M = size (F,1).
FORTRAN 90 Interface
Generic: CALL GCIN (P, X, F [, …])
Specific: The specific interface names are S_GCIN and D_GCIN.
FORTRAN 77 Interface
Single: CALL GCIN (P, IOPT, M, X, F)
Double: The double precision function name is DGCIN.
Description
Function GCIN evaluates the inverse of a continuous distribution function, given ordinates of the probability density function. The range of the distribution must be specified in X. For distributions with infinite ranges, endpoints must be chosen so that most of the probability content is included.
The function GCIN first fits a curve to the points given in X and F with either a piecewise linear interpolant or a C1 cubic spline interpolant based on a method by Akima (1970). Function GCIN then determines the area, A, under the curve. (If the distribution were of finite range and if the fit were exact, this area would be 1.0.) It next finds the maximum abscissa up to which the area is less than AP and the minimum abscissa up to which the area is greater than AP. The routine then interpolates for the point corresponding to AP. Because of the scaling by A, it is not assumed that the integral of the density defined by X and F is 1.0.
For most distributions, it is likely that better approximations to the distribution function are obtained when IOPT equals 3 or 4, that is, when a cubic spline is used to approximate the function. It is also likely that better approximations can be obtained when the abscissas are chosen more densely over regions where the density and its derivatives (when they exist) are varying greatly.
Comments
Workspace may be explicitly provided, if desired, by the use of G3IN/DG3IN. The reference is
G3IN(P, IOPT, M, X, F, WK, IWK)
The arguments in addition to those of GCIN are:
WK — Work vector of length 5 * M if IOPT = 3, and of length 4 * M if IOPT = 4.
IWK — Work vector of length M.
Example
In this example, we find the 90‑th percentage point for a beta random variable with parameters 12 and 12. The probability density function of a beta random variable with parameters p and q is
where
Γ(
⋅) is the gamma function. The density is equal to 0 outside the interval [0, 1]. With
p =
q, this is a symmetric distribution. Knowing that the probability density of this distribution is very peaked in the vicinity of 0.5, we could perhaps get a better fit by using unequally spaced abscissas, but we will keep it simple and use 300 equally spaced points. Note that this is the same example that is used in the description of routine
BETIN. The result from
BETIN would be expected to be more accurate than that from
GCIN since
BETIN is designed specifically for this distribution.
USE GCIN_INT
USE UMACH_INT
USE BETA_INT
IMPLICIT NONE
INTEGER M
PARAMETER (M=300)
!
INTEGER I, IOPT, NOUT
REAL C, F(M), H, P, PIN, PIN1, QIN, QIN1, &
X(2), X0, XI
!
CALL UMACH (2, NOUT)
P = 0.9
IOPT = 3
! Initializations for a beta(12,12)
! distribution.
PIN = 12.0
QIN = 12.0
PIN1 = PIN - 1.0
QIN1 = QIN - 1.0
C = 1.0/BETA(PIN,QIN)
XI = 0.0
H = 1.0/(M-1.0)
X(1) = XI
F(1) = 0.0
XI = XI + H
! Compute ordinates of the probability
! density function.
DO 10 I=2, M - 1
F(I) = C*XI**PIN1*(1.0-XI)**QIN1
XI = XI + H
10 CONTINUE
X(2) = 1.0
F(M) = 0.0
X0 = GCIN(P,X,F, IOPT=IOPT)
WRITE (NOUT,99999) X0
99999 FORMAT (' X is less than ', F6.4, ' with probability 0.9.')
END
Output
X is less than 0.6304 with probability 0.9.
