R/qgrubbs.R

In outliers: Tests for Outliers

"qgrubbs" <-
function(p,n,type=10,rev=FALSE)
{

if (type == 10) {

if (!rev) { 
return(((n-1)/sqrt(n))*sqrt(qt((1-p)/n,n-2)^2/(n-2+qt((1-p)/n,n-2)^2)))
}
else{
s <- (p^2*n*(2-n))/(p^2*n-(n-1)^2)
t <- sqrt(s)
if (is.nan(t)) {
        res <- 0
}
else
{
res <- n*(1-pt(t,n-2))
res[res>1] <- 1
}
return(1-res)
}
}

else if (type == 11) {

if (!rev) {
return(sqrt((2*(n-1)*qt((1-p)/(n*(n-1)),n-2)^2)/(n-2+qt((1-p)/(n*(n-1)),n-2)^2)))
}
else{

q <- p;

p <- vector();

for (i in 1:length(q)) {

if (q[i] > qgrubbs(0.9999,n,type=11)) { pp <- 1 }
else
if (q[i] < qgrubbs(2e-16,n,type=11)) { pp <- 0}

else {
        f <- function(x,q,n) { qgrubbs(x,n,type=11)-q }
        pp <- uniroot(f,c(0.001,0.9999),q=q[i],n=n)$root

}
p <- c(p,pp)

}

return(p)

}

}

else {

if (n > 30) stop("n must be in range 3-30");

pp <- c(0.01,0.025,0.05,0.1,0.15,0.2,0.4,0.6,0.8,0.9,0.95,0.975,0.99);

gtwo <- c(
NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
0.00001,0.0002,0.0008,0.0031,0.007,0.013,0.055,0.138,0.283,0.399,0.482,0.54,0.589,
0.0035,0.009,0.0183,0.0376,0.058,0.078,0.169,0.276,0.41,0.502,0.571,0.63,0.689,
0.0186,0.0349,0.0565,0.0921,0.124,0.153,0.257,0.361,0.478,0.562,0.626,0.674,0.724,
0.044,0.0708,0.102,0.1479,0.186,0.217,0.325,0.423,0.53,0.605,0.659,0.701,0.743,
0.075,0.1101,0.1478,0.1994,0.238,0.271,0.376,0.468,0.568,0.635,0.684,0.722,0.763,
0.1082,0.1492,0.1909,0.2454,0.285,0.319,0.42,0.507,0.598,0.658,0.703,0.738,0.776,
0.1415,0.1865,0.2305,0.2863,0.326,0.358,0.455,0.537,0.622,0.678,0.721,0.753,0.787,
0.1736,0.2212,0.2666,0.3226,0.363,0.394,0.487,0.564,0.624,0.695,0.734,0.765,0.797,
0.2044,0.2536,0.2996,0.3552,0.393,0.424,0.514,0.585,0.66,0.709,0.745,0.773,0.802,
0.2333,0.2836,0.3295,0.3843,0.421,0.451,0.537,0.605,0.676,0.721,0.755,0.782,0.809,
0.2605,0.3112,0.3568,0.4106,0.447,0.477,0.558,0.622,0.689,0.733,0.765,0.789,0.817,
0.2859,0.3367,0.3818,0.4345,0.469,0.497,0.576,0.639,0.701,0.742,0.773,0.797,0.822,
0.3098,0.3603,0.4048,0.4562,0.491,0.518,0.593,0.653,0.713,0.751,0.78,0.803,0.826,
0.3321,0.3822,0.4259,0.4761,0.511,0.536,0.609,0.666,0.723,0.76,0.787,0.808,0.831,
0.353,0.4025,0.4455,0.4944,0.526,0.552,0.622,0.676,0.732,0.767,0.793,0.814,0.835,
0.3725,0.4214,0.4636,0.5113,0.543,0.567,0.635,0.688,0.74,0.774,0.799,0.818,0.838,
0.3909,0.4391,0.4804,0.5269,0.559,0.582,0.647,0.697,0.748,0.781,0.805,0.823,0.843,
0.408,0.457,0.496,0.542,0.571,0.594,0.657,0.706,0.755,0.786,0.81,0.828,0.847,
0.425,0.474,0.512,0.556,0.584,0.606,0.668,0.715,0.762,0.792,0.815,0.834,0.85,
0.442,0.486,0.524,0.568,0.596,0.618,0.677,0.723,0.769,0.797,0.819,0.836,0.853,
0.453,0.5,0.538,0.581,0.608,0.628,0.686,0.73,0.774,0.802,0.823,0.84,0.857,
0.466,0.511,0.547,0.589,0.616,0.637,0.693,0.736,0.779,0.807,0.827,0.843,0.86,
0.482,0.525,0.561,0.601,0.627,0.647,0.701,0.743,0.784,0.811,0.83,0.845,0.861,
0.492,0.536,0.572,0.611,0.636,0.655,0.709,0.749,0.789,0.815,0.834,0.849,0.864,
0.505,0.548,0.583,0.621,0.646,0.664,0.716,0.755,0.794,0.819,0.837,0.851,0.866,
0.516,0.558,0.592,0.629,0.654,0.672,0.722,0.76,0.798,0.822,0.84,0.854,0.869,
0.528,0.568,0.602,0.638,0.661,0.679,0.728,0.765,0.802,0.826,0.842,0.856,0.87);

dim(gtwo) <- c(13,30);

if (!rev) res <- qtable(p,pp,gtwo[,n])
else res <- qtable(p,gtwo[,n],pp)

res[res < 0] <- 0;
res[res > 1] <- 1;


return(res)

}
}