EMS {MGBT} | R Documentation |
Expected values of M and S
Description
Compute expected values of M
and S
given q_\mathrm{min}
and define the quantity
z_r = \Phi^{(1)}(q_\mathrm{min})\mbox{,}
where \Phi^{(1)}(\cdot)
is the inverse of the standard normal distribution. As result, q_\mathrm{min}
is itself a probability because it is an argument to the qnorm()
function. The expected value M
is defined as
M = \Psi(z_r, 1)\mbox{,}
where \Psi(a,b)
is the gtmoms
function. The S
requires the conditional moments of the Chi-square (CondMomsChi2
) defined as the two value vector \mbox{}_2S
that provides the values \alpha = \mbox{}_2S_1^2 / \mbox{}_2S_2
and \beta = \mbox{}_2S_2 / \mbox{}_2S_1
. The S
is then defined by
S = \sqrt{\beta}\biggl(\frac{\Gamma(\alpha+0.5)}{\Gamma(\alpha)}\biggr)\mbox{.}
Usage
EMS(n, r, qmin)
Arguments
n |
The number of observations; |
r |
The number of truncated observations? (confirm); and |
qmin |
A nonexceedance probability threshold for |
Value
The expected values of M
and S
in the form of an R vector
.
Note
TAC sources call on the explicit first index of M
as literally “Em[1]
” for the returned vector, which seems unnecessary. This is a potential weak point in design because the gtmoms
function is naturally vectorized and could potentially produce a vector of M
values. For the implementation here, only the first value in qmin
is used and a warning otherwise issued. Such coding prevents the return value from EMS
accidentally acquiring a length greater than two. For at least samples of size n=2
, overranging in a call to lgamma(alpha)
happens for alpha=0
. A suppressWarnings()
is wrapped around the applicable line. The resulting NaN
cascades up the chain, which will end up inside peta
, but therein SigmaMp
is not finite and a p-value of unity is returned.
Author(s)
W.H. Asquith consulting T.A. Cohn sources
Source
LowOutliers_jfe(R).txt
, LowOutliers_wha(R).txt
, P3_089(R).txt
—Named EMS
References
Cohn, T.A., 2013–2016, Personal communication of original R source code: U.S. Geological Survey, Reston, Va.
See Also
CondMomsChi2
, EMS
, VMS
, V
, gtmoms
Examples
EMS(58,2,.5)
#[1] 0.7978846 0.5989138
# Monte Carlo experiment to test EMS and VMS functions
"test_EMS" <- function(nrep=1000, n=100, r=0, qr=0.2, ss=1) { # TAC named function
set.seed(ss)
Moms <- replicate(n=nrep, {
x <- qnorm(runif(n-r,min=qr,max=1));
c(mean(x),var(x))}); xsi <- qnorm(qr);
list(
MeanMS_obs = c(mean(Moms[1,]), mean(sqrt(Moms[2,])), mean(Moms[2,])),
EMS = c(EMS(n,r,qr), gtmoms(xsi,2) - gtmoms(xsi,1)^2),
CovMS2_obs = cov(t(Moms)),
VMS2 = V(n,r,qr),
VMS_obs = array(c(var( Moms[1,]),
rep(cov( Moms[1,], sqrt(Moms[2,])),2),
var(sqrt(Moms[2,]))), dim=c(2,2)),
VMS = VMS(n,r,qr) )
}
test_EMS()