expect.max.ostat: Compute the Expectation of a Maximum (or Minimum and others) Order Statistic

Description

The maximum (or minimum) expectation of an order statistic can be directly used for L-moment computation through either of the following two equations (Hosking, 2006) as dictated by using the maximum ($\mathrm{E}[X_{k:k}]$, expect.max.ostat) or minimum ($\mathrm{E}[X_{1:k}]$, expect.min.ostat): $$ \lambda_r = (-1)^{r-1} \sum_{k=1}^r (-1)^{r-k}k^{-1}{r-1 \choose k-1}{r+k-2 \choose k-1}\mathrm{E}[X_{1:k}]\mbox{,} $$ and $$ \lambda_r = \sum_{k=1}^r (-1)^{r-k}k^{-1}{r-1 \choose k-1}{r+k-2 \choose k-1}\mathrm{E}[X_{k:k}]\mbox{.} $$

In terms of the quantile function qlmomco, the expectation of an order statistic (Asquith, 2011, p. 49) is $$ \mathrm{E}[X_{j:n}] = n {n-1 \choose j - 1}\int^1_0 \! x(F)\times F^{j-1} \times (1-F)^{n-j}\; \mathrm{d}F\mbox{,} $$ where $x(F)$ is the quantile function, $F$ is nonexceedance probability, $n$ is sample size, and $j$ is the $j$th order statistic.

In terms of the probability density function (PDF) dlmomco and cumulative density function (CDF) plmomco, the expectation of an order statistic (Asquith, 2011, p. 50) is $$ \mathrm{E}[X_{j:n}] = \frac{1}{\mathrm{B}(j,n-j+1)}\int_{-\infty}^{\infty} [F(x)]^{j-1}[1-F(x)]^{n-j} x\, f(x)\;\mathrm{d} x\mbox{,} $$ where $F(x)$ is the CDF, $f(x)$ is the PDF, and $\mathrm{B}(j,n-j+1)$ is the complete Beta function, which in R is beta with the same argument order as shown above.

Usage

expect.max.ostat(n, para=NULL, cdf=NULL, pdf=NULL, qua=NULL,
                 j=NULL, lower=-Inf, upper=Inf, aslist=FALSE, ...)

Arguments

The sample size.

para

A distribution parameter list from a function such as vec2par or lmom2par.

cdf

cumulative distribution function of the distribution.

pdf

probability density function of the distribution.

qua

quantile function of the distribution. If this is defined, then cdf and pdf are ignored.

The $j$th value of the order statistic, which defaults to n=j (the maximum order statistic) if j=NULL.

lower

The lower limit for integration.

upper

The upper limit for integration.

aslist

A logically triggering whether an R list is returned instead of just the expection.

...

Additional arguments to pass to the three distribution functions.

Value

The expectation of the maximum order statistic, unless $j$ is specified and then the expectation of that order statistic is returned. This similarly holds if the expect.min.ostat function is used except “maximum” becomes the “minimum”.

Alternatively, an R list is returned.

type

The type of approach used: “bypdfcdf” means the PDF and CDF of the distribution were used, and alternatively “byqua” means that the quantile function was used.

value

See previous discussion of value.

abs.error

Estimate of the modulus of the absolute error from R function integrate.

subdivisions

The number of subintervals produced in the subdivision process from R function integrate.

message

“OK” or a character string giving the error message.

Details

If qua != NULL, then the first order-statistic expectation equation above is used, and any function that might have been set in cdf and pdf is ignored. If the limits are infinite (default), then the limits of the integration will be set to $F\!\downarrow = 0$ and $F\!\uparrow = 1$. The user can replace these by setting the limits to something “near” zero and(or) “near” 1. Please consult the Note below concerning more information about the limits of integration.

If qua == NULL, then the second order-statistic expectation equation above is used and cdf and pdf must be set. The default $\pm\infty$ limits are used unless the user knows otherwise for the distribution or through supervision provides their meaning of small and large.

This function requires the user to provide either the qua or the cdf and pdf functions, which is somewhat divergent from the typical flow of logic of lmomco. This has been done so that expect.max.ostat can be used readily for experimental distribution functions. It is suggested that the parameter object be left in the lmomco style (see vec2par) even if the user is providing their own distribution functions.

Last comments: This function is built around the idea that either (1) the cdf and pdf ensemble or (2) qua exist in some clean analytical form and therefore the qua=NULL is the trigger on which order statistic expectation integral is used. This precludes an attempt to compute the support of the distribution internally, and thus providing possibly superior (more refined) lower and upper limits. Here is a suggested re-implementation using the support of the Generalized Extreme Value distribution:

para <- vec2par(c(100, 23, -0.5), type="gev")
lo <- quagev(0, para) # The value 54
hi <- quagev(1, para) # Infinity
E22 <- expect.max.ostat(2, para=para,cdf=cdfgev,pdf=pdfgev,
                           lower=lo, upper=hi)
E21 <- expect.min.ostat(2, para=para,cdf=cdfgev,pdf=pdfgev,
                           lower=lo, upper=hi)
L2 <- (E22 - E21)/2 # definition of L-scale
cat("L-scale: ",L2,"(integration)",
    lmomgev(para)$lambdas[2], "(theory)\n")
# The results show 33.77202 as L-scale.

The design intent makes it possible for some arbitrary and(or) new quantile function with difficult cdf and pdf expressions (or numerical approximations) to not be needed as the L-moments are explored. Contrarily, perhaps some new pdf exists and simple integration of it is made to get the cdf but the qua would need more elaborate numerics to invert the cdf. The user could then still explore the L-moments with supervision on the integration limits or foreknowledge of the support of the distribution.

References

Asquith, W.H., 2011, Distributional analysis with L-moment statistics using the R environment for statistical computing: Createspace Independent Publishing Platform, ISBN 978--146350841--8.

Gilchrist, W.G., 2000, Statistical modelling with quantile functions: Chapman and Hall/CRC, Boca Raton.

Hosking, J.R.M., 2006, On the characterization of distributions by their L-moments: Journal of Statistical Planning and Inference, v. 136, no. 1, pp. 193--198.

Examples

Run this code

# NOT RUN {
para <- vec2par(c(10,100), type="nor")
n <- 12
# The three outputted values from should be similar:
# (1) theoretical, (2) theoretical, and (3) simulation
expect.max.ostat(n, para=para, cdf=cdfnor, pdf=pdfnor)
expect.max.ostat(n, para=para, qua=quanor)
mean(sapply(1:1000, function(x) { max(rlmomco(n,para))}))

eostat(8, j=5, qua=quagum, para=vec2par(c(1670,1000), type="gum"))
# }
# NOT RUN {
para <- vec2par(c(1280, 800), type="nor")
expect.max.ostat(10, j=9, para, qua=quanor)
[1] 2081.086      # SUCCESS ---------------------------
expect.max.ostat(10, j=9, para, pdf=pdfnor, cdf=cdfnor,
                                lower=-1E3, upper=1E6)
[1] 1.662701e-06  # FAILURE ---------------------------
expect.max.ostat(10, j=9, para, pdf=pdfnor, cdf=cdfnor,
                                lower=-1E3, upper=1E5)
[1] 2081.086      # SUCCESS ---------------------------
# }

Run the code above in your browser using DataLab