Evaluates the inverse of a general continuous cumulative distribution function given ordinates of the density.
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.
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).
Generic: CALL GCIN (P, X, F [,…])
Specific: The specific interface names are S_GCIN and D_GCIN.
Single: CALL GCIN (P, IOPT, M, X, F)
Double: The double precision function name is DGCIN.
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.
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.
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
X is less than 0.6304 with probability 0.9.
PHONE: 713.784.3131 FAX:713.781.9260 |