R/CEPs.R
In dwoll/shotGroups: Analyze Shot Group Data
Defines functions
CEPValstar
CEPRAND
CEPEthridge
CEPRMSE
CEPIgnani
CEPKrempasky
CEPZhang
CEPRayleigh
CEPGrubbsLiu
CEPGrubbsPearson
CEPGrubbsPatnaik
CEPCorrNormal
#####-----------------------------------------------------------------------
## CEP estimate based on correlated bivariate normal distribution - mvnEll.R, hoyt.R
CEPCorrNormal <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
CorrNormal <- if(accuracy) { # POA != POI
## quantile from offset circle probability - most general case
qmvnEll(CEPlevel, mu=numeric(p), sigma=sigma, e=diag(p), x0=ctr)
} else { # POA == POI
if(p == 2L) { # 2D -> exact Hoyt distribution
HP <- getHoytParam(sigma)
qHoyt(CEPlevel, qpar=HP$q, omega=HP$omega)
} else { # 1D/3D case
qmvnEll(CEPlevel, mu=numeric(p), sigma=sigma, e=diag(p), x0=numeric(p))
}
}
setNames(CorrNormal, CEPlevel)
}
####-----------------------------------------------------------------------
## Grubbs estimates (Grubbs, 1964, p54, p55-56) - grubbs.R
## Grubbs-Patnaik CEP estimate based on Patnaik two-moment central chi^2 approximation
CEPGrubbsPatnaik <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE) {
GPP <- getGrubbsParam(sigma, ctr=ctr, accuracy=accuracy)
GrubbsPatnaik <- qChisqGrubbs(CEPlevel, m=GPP$m, v=GPP$v, n=GPP$n, type="Patnaik")
setNames(GrubbsPatnaik, CEPlevel)
}
## Grubbs-Pearson CEP estimate based on Pearson three-moment central chi^2 approximation
CEPGrubbsPearson <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE) {
GPP <- getGrubbsParam(sigma, ctr=ctr, accuracy=accuracy)
GrubbsPearson <- qChisqGrubbs(CEPlevel, m=GPP$m, v=GPP$v, nPrime=GPP$nPrime, type="Pearson")
setNames(GrubbsPearson, CEPlevel)
}
## Grubbs-Liu CEP estimate based on four-moment non-central chi^2
## approximation (Liu, Tang & Zhang, 2009)
CEPGrubbsLiu <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE) {
GPP <- getGrubbsParam(sigma, ctr=ctr, accuracy=accuracy)
GrubbsLiu <- qChisqGrubbs(CEPlevel, m=GPP$m, v=GPP$v, muX=GPP$muX,
varX=GPP$varX, l=GPP$l, delta=GPP$delta, type="Liu")
setNames(GrubbsLiu, CEPlevel)
}
#####-----------------------------------------------------------------------
## Rayleigh CEP estimate from Singh, 1992 - rayleigh.R, rice.R, maxwell.R
CEPRayleigh <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE, doRob=FALSE, xy) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
Rayleigh <- if(p == 1L) { # 1D
if(accuracy) { # POA != POI -> non-central F distribution
RiceParam <- getRiceParam(xy, doRob=doRob)
N <- length(xy)
NU <- RiceParam$nu
SD <- setNames(RiceParam$sigma["sigma"], NULL)
NCP <- 1*(NU/SD)^2 # N*f^2, N=1 because interested in individual shot, not mean
SD * sqrt(qf(CEPlevel, df1=1, df2=N-1, ncp=NCP))
} else { # POA = POI -> half normal / central F distribution
N <- length(xy)
SD <- getRayParam(xy, doRob=doRob)$sigma["sigma"]
SD * sqrt(qchisq(CEPlevel, df=p))
}
} else if(p == 2L) { # 2D
if(accuracy) { # POA != POI -> Rice distribution
RiceParam <- getRiceParam(xy, doRob=doRob)
qRice(CEPlevel, nu=RiceParam$nu, sigma=RiceParam$sigma["sigma"])
} else { # POA = POI -> Rayleigh
RayParam <- getRayParam(xy, doRob=doRob)
qRayleigh(CEPlevel, scale=RayParam$sigma["sigma"])
}
} else if(p == 3L) { # 3D
MaxParam <- getRayParam(xy, doRob=doRob)
if(accuracy) { # POA != POI -> offset sphere probability
## circular covariance matrix with estimated M-B param sigma
sigMat <- diag(rep(MaxParam$sigma["sigma"]^2, p))
qmvnEll(CEPlevel, sigma=sigMat, mu=numeric(p), x0=ctr, e=diag(p))
} else { # POA = POI -> Maxwell-Boltzmann distribution
qMaxwell(CEPlevel, sigma=MaxParam$sigma["sigma"])
}
} else {
warning("Rayleigh CEP is only available for 1D/2D/3D-data")
rep(NA_real_, length(CEPlevel))
}
setNames(Rayleigh, CEPlevel)
}
#####-----------------------------------------------------------------------
## Zhang estimate from Zhang & An (2012) p570 eq. 7
CEPZhang <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE,
isotropic=FALSE, doRob=FALSE, xy) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
## make sure eigenvalues >= 0 when very small
ev <- eigen(sigma)$values # eigenvalues
lambda <- ev*sign(ev)
if(!accuracy) {
ctr <- c(0, 0)
}
Zhang50 <- if((p == 2L) && isotropic) { # assume equal variances
sigma <- getRayParam(xy=xy, level=CEPlevel, doRob=doRob)[["sigma"]]["sigma"]
r0 <- sqrt(sum((ctr/sigma)^2)) # p 570
print(r0)
if(r0 <= sqrt(5)) {
A0 <- 1.17712063
A1 <- 0.01635468
A2 <- 0.21387229
A3 <- 0.13972335
A4 <- -0.08426949
A5 <- 0.01284512
sigma*(A0 + A1*r0 + A2*r0^2 + A3*r0^3 + A4*r0^4 + A5*r0^5)
} else {
if(r0 <= 3) {
d <- 4
B0 <- 2.2458024
B1 <- 0.2200099
B2 <- -0.01003873
B3 <- 0.00078179
B4 <- -0.00005397
B5 <- 0.00000204
} else if(r0 <= 4) {
d <- 9
B0 <- 3.165246
B1 <- 0.15761986
B2 <- -0.00387106
B3 <- 0.00017928
B4 <- -0.00000822
B5 <- 0.00000022
} else if(r0 <= 5) {
d <- 16
B0 <- 4.1243779
B1 <- 0.12114747
B2 <- -0.00177091
B3 <- 0.00005016
B4 <- -0.0000015
B5 <- 0.00000003
} else if(r0 <= sqrt(45)) {
d <- 25
B0 <- 5.099676
B1 <- 0.09800860
B2 <- -0.00093530
B3 <- 0.00001679
B4 <- -0.00000029
B5 <- 0
} else {
warning("Isotropic case only defined up to r0 <= sqrt(45)")
d <- NA_real_
B0 <- NA_real_
B1 <- NA_real_
B2 <- NA_real_
B3 <- NA_real_
B4 <- NA_real_
B5 <- NA_real_
}
sigma * (B0 + B1*(r0^2 - d) +
B2*(r0^2 - d)^2 +
B3*(r0^2 - d)^3 +
B4*(r0^2 - d)^4 +
B5*(r0^2 - d)^5)
}
} else if((p == 2L)) { # do not assume equal variances
r1 <- sqrt(sum(ctr^2/lambda)) # p568
## coefficients for eq. 7 (personal communication)
Bijk <- c("i j k B
0 0 0 -1.50800
0 0 1 37.06914
0 0 2 -261.85730
0 0 3 992.70823
0 0 4 -2127.35633
0 0 5 2432.20862
0 0 6 -866.72077
0 0 7 -1043.20376
0 0 8 1007.57426
0 0 9 201.09441
0 0 10 -547.24802
0 0 11 178.38609
0 1 0 30.73407
0 1 1 -559.72454
0 1 2 4340.04504
0 1 3 -18573.68822
0 1 4 48237.39491
0 1 5 -78643.98158
0 1 6 80365.58824
0 1 7 -51299.58759
0 1 8 24063.13327
0 1 9 -14082.77353
0 1 10 8340.22172
0 1 11 -2217.07894
0 2 0 -140.53635
0 2 1 2615.24869
0 2 2 -21136.32531
0 2 3 95616.46894
0 2 4 -266698.90970
0 2 5 478107.64088
0 2 6 -557690.04739
0 2 7 425021.90498
0 2 8 -221989.27415
0 2 9 95629.50576
0 2 10 -38007.28570
0 2 11 8670.58029
0 3 0 307.16842
0 3 1 -5592.48090
0 3 2 45454.19525
0 3 3 -210124.94755
0 3 4 605819.65479
0 3 5 -1134958.18022
0 3 6 1399681.23997
0 3 7 -1135207.32583
0 3 8 611593.91350
0 3 9 -236561.30418
0 3 10 74259.96471
0 3 11 -14670.14344
0 4 0 -338.98709
0 4 1 5886.07371
0 4 2 -46868.89957
0 4 3 216011.08369
0 4 4 -628036.11139
0 4 5 1196820.36257
0 4 6 -1511973.91364
0 4 7 1259646.73040
0 4 8 -686657.86707
0 4 9 251291.10861
0 4 10 -67419.94690
0 4 11 11638.94476
0 5 0 150.64942
0 5 1 -2481.96624
0 5 2 19094.1026
0 5 3 -86299.71627
0 5 4 248634.95180
0 5 5 -473172.89850
0 5 6 600573.05422
0 5 7 -504295.29927
0 5 8 275477.36039
0 5 9 -97818.1152801
0 5 10 23834.43367
0 5 11 -3696.11403
1 0 0 7.30917
1 0 1 -121.95931
1 0 2 840.42710
1 0 3 -3062.11230
1 0 4 6186.79357
1 0 5 -6054.59489
1 0 6 -208.97884
1 0 7 6535.71511
1 0 8 -5133.24462
1 0 9 -260.95477
1 0 10 1917.57177
1 0 11 -645.78090
1 1 0 -110.31830
1 1 1 2017.70282
1 1 2 -15683.19635
1 1 3 67488.60741
1 1 4 -176701.28435
1 1 5 290580.88672
1 1 6 -298955.31158
1 1 7 191453.37431
1 1 8 -90293.34635
1 1 9 54007.00700
1 1 10 -32486.81832
1 1 11 8681.02815
1 2 0 498.12492
1 2 1 -9427.49801
1 2 2 77155.89222
1 2 3 -354370.79904
1 2 4 10006477.27304
1 2 5 -1840218.83052
1 2 6 2190104.05741
1 2 7 -1701102.53064
1 2 8 899284.68715
1 2 9 -383894.11296
1 2 10 149222.77548
1 2 11 -33723.07562
1 3 0 -1035.26688
1 3 1 19430.28736
1 3 2 -161868.38075
1 3 3 767936.17634
1 3 4 -2276698.65746
1 3 5 4389361.94030
1 3 6 -5567199.26107
1 3 7 4629872.46270
1 3 8 -2529969.48968
1 3 9 961497.85815
1 3 10 -285581.48386
1 3 11 54243.94576
1 4 0 1101.18523
1 4 1 -19807.40970
1 4 2 162886.78484
1 4 3 -776205.30786
1 4 4 2336828.66380
1 4 5 -4612153.15252
1 4 6 6025709.64684
1 4 7 -5168691.92531
1 4 8 2862608.29402
1 4 9 -1024135.74071
1 4 10 251524.75518
1 4 11 -39658.01666
1 5 0 -487.75337
1 5 1 8290.53777
1 5 2 -65870.58342
1 5 3 308197.20491
1 5 4 -921057.63790
1 5 5 1819395.02951
1 5 6 -2393829.33582
1 5 7 2073889.98097
1 5 8 -1152239.51773
1 5 9 398819.42046
1 5 10 -86623.36549
1 5 11 11513.64407
2 0 0 -8.39680
2 0 1 137.59442
2 0 2 -926.41844
2 0 3 3265.94858
2 0 4 -6100.14231
2 0 5 4278.05349
2 0 6 4933.42238
2 0 7 -12641.00423
2 0 8 8961.99007
2 0 9 -256.37747
2 0 10 -2517.52899
2 0 11 872.69812
2 1 0 151.85852
2 1 1 -2764.43629
2 1 2 21418.68153
2 1 3 -92079.73132
2 1 4 240640.93680
2 1 5 -392722.30429
2 1 6 395837.40643
2 1 7 -243259.82654
2 1 8 112039.90079
2 1 9 -74462.27711
2 1 10 48378.10091
2 1 11 -13174.53964
2 2 0 -696.05269
2 2 1 13235.79923
2 2 2 -108719.94280
2 2 3 502418.93939
2 2 4 -1437223.47953
2 2 5 2643175.11851
2 2 6 -3154239.38546
2 2 7 2448763.15222
2 2 8 -1297826.56619
2 2 9 566806.86704
2 2 10 -228178.84490
2 2 11 52471.09248
2 3 0 1400.85858
2 3 1 -26862.13543
2 3 2 227395.73358
2 3 3 -1096233.13091
2 3 4 3302421.47593
2 3 5 -6460706.72431
2 3 6 8294685.64827
2 3 7 -6962294.81091
2 3 8 3829716.07859
2 3 9 -1463576.88726
2 3 10 437814.27491
2 3 11 -83738.96835
2 4 0 -1432.43625
2 4 1 26735.89822
2 4 2 -226271.58166
2 4 3 1106145.22484
2 4 4 -3409574.42966
2 4 5 6873381.61298
2 4 6 -9144047.77589
2 4 7 7956459.96580
2 4 8 -4444585.05705
2 4 9 1586114.68543
2 4 10 -381893.75700
2 4 11 58951.23793
2 5 0 622.91812
2 5 1 -11035.33184
2 5 2 90838.15831
2 5 3 -438500.09965
2 5 4 1348016.12322
2 5 5 -2731098.77184
2 5 6 3673259.51515
2 5 7 -3238879.5041
2 5 8 1817714.84110
2 5 9 -624450.71184
2 5 10 129495.99917
2 5 11 -15978.20880
3 0 0 4.80967
3 0 1 -77.99167
3 0 2 516.96171
3 0 3 -1760.05286
3 0 4 2978.06661
3 0 5 -938.44505
3 0 6 -5613.89495
3 0 7 10344.99174
3 0 8 -6992.49989
3 0 9 -507.98776
3 0 10 1593.73407
3 0 11 -563.14246
3 1 0 -106.25020
3 1 1 1916.56093
3 1 2 -14700.28712
3 1 3 762447.47246
3 1 4 -160293.41157
3 1 5 253462.21660
3 1 6 -240535.58640
3 1 7 131669.95187
3 1 8 -55767.81377
3 1 9 47198.42290
3 1 10 -35116.21231
3 1 11 9820.81246
3 2 0 511.67374
3 2 1 -9673.11429
3 2 2 78791.45853
3 2 3 -360766.93776
3 2 4 1020321.00332
3 2 5 -1846866.00449
3 2 6 2154711.78200
3 2 7 -1626738.07815
3 2 8 852404.97838
3 2 9 -395825.66688
3 2 10 175015.53634
3 2 11 -41872.09886
3 3 0 -1027.83063
3 3 1 19916.0826
3 3 2 -168649.08473
3 3 3 810228.93041
3 3 4 -2426699.32685
3 3 5 4707457.17127
3 3 6 -5975393.05967
3 3 7 4953340.57762
3 3 8 -2714680.08066
3 3 9 1072562.25276
3 3 10 -347415.58697
3 3 11 70336.92513
3 4 0 1011.50846
3 4 1 -19520.35600
3 4 2 167740.84261
3 4 3 -825219.08597
3 4 4 2547613.20509
3 4 5 -5126865.40135
3 4 6 6791764.92505
3 4 7 -5879408.35900
3 4 8 8282864.52405
3 4 9 -1196814.35705
3 4 10 308131.74625
3 4 11 -51281.82176
3 5 0 -419.17227
3 5 1 7799.23359
3 5 2 -66151.92153
3 5 3 325082.85879
3 5 4 -1009967.60753
3 5 5 2058208.63394
3 5 6 -2775472.43675
3 5 7 2449243.21250
3 5 8 -1377607.57999
3 5 9 479600.83760
3 5 10 -104225.66822
3 5 11 13904.43044
4 0 0 -1.41012
4 0 1 22.70521
4 0 2 -148.46582
4 0 3 489.74223
4 0 4 -747.89192
4 0 5 -89.38056
4 0 6 2297.32756
4 0 7 -3725.25391
4 0 8 2452.57493
4 0 9 -251.29178
4 0 10 -464.39654
4 0 11 165.46501
4 1 0 36.91482
4 1 1 -657.65046
4 1 2 4965.58653
4 1 3 -20654.80936
4 1 4 51363.19899
4 1 5 -76954.43695
4 1 6 65622.22571
4 1 7 -27839.71942
4 1 8 9047.66734
4 1 9 -13506.39268
4 1 10 12032.45580
4 1 11 -3452.84254
4 2 0 -189.23560
4 2 1 3534.90907
4 2 2 -28340.16312
4 2 3 127301.53392
4 2 4 -351465.57188
4 2 5 616163.17964
4 2 6 -688181.61871
4 2 7 492046.07078
4 2 8 -251414.70516
4 2 9 129307.34993
4 2 10 -65040.55204
4 2 11 16270.97692
4 3 0 389.49639
4 3 1 -7517.54072
4 3 2 62883.56972
4 3 3 -297233.52178
4 3 4 872933.18776
4 3 5 -1653991.63581
4 3 6 2041991.02768
4 3 7 -1644559.84876
4 3 8 891362.10361
4 3 9 -372834.21958
4 3 10 136205.14245
4 3 11 -29615.35099
4 4 0 -375.58749
4 4 1 7349.92108
4 4 2 -63036.60176
4 4 3 307356.35957
4 4 4 -936607.41424
4 4 5 1854763.70135
4 4 6 -2412369.0135
4 4 7 2051411.99056
4 4 8 -1138662.20737
4 4 9 431751.74719
4 4 10 -124649.03235
4 4 11 23056.76527
4 5 0 146.79700
4 5 1 -2827.54869
4 5 2 24329.22818
4 5 3 -120032.35962
4 5 4 372079.75852
4 5 5 -753454.52592
4 5 6 1006974.01091
4 5 7 -880553.11272
4 5 8 494280.92111
4 5 9 -176872.67948
4 5 10 42410.61903
4 5 11 -6478.29063
5 0 0 0.16415
5 0 1 -2.62846
5 0 2 17.01287
5 0 3 -54.90894
5 0 4 78.17490
5 0 5 32.16035
5 0 6 -297.35896
5 0 7 454.60877
5 0 8 -293.38420
5 0 9 -37.40186
5 0 10 44.46108
5 0 11 -15.65335
5 1 0 -4.95117
5 1 1 86.91975
5 1 2 -643.82509
5 1 3 2608.72592
5 1 4 -6236.22491
5 1 5 8724.68903
5 1 6 -6379.03883
5 1 7 1475.79282
5 1 8 45.93268
5 1 9 1381.76407
5 1 10 -1495.997
5 1 11 435.75596
5 2 0 27.39345
5 2 1 -503.58840
5 2 2 3956.27375
5 2 3 -17332.34425
5 2 4 46350.38848
5 2 5 -77842.78318
5 2 6 81803.47877
5 2 7 -53864.40054
5 2 8 26400.02294
5 2 9 -15729.73024
5 2 10 9100.16729
5 2 11 -2363.21625
5 3 0 -58.88322
5 3 1 1121.35889
5 3 2 -9198.10897
5 3 3 42452.68426
5 3 4 -121199.98901
5 3 5 221998.94268
5 3 6 -263183.30573
5 3 7 202902.59108
5 3 8 -108017.59066
5 3 9 48909.67004
5 3 10 -20479.54143
5 3 11 4749.49844
5 4 0 56.60322
5 4 1 -1105.18410
5 4 2 961.78170
5 4 3 -44835.28439
5 4 4 133645.11524
5 4 5 -257933.62666
5 4 6 325955.13839
5 4 7 -269490.82354
5 4 8 148047.67616
5 4 9 -59284.13086
5 4 10 19575.20456
5 4 11 -3990.41775
5 5 0 -20.77692
5 5 1 407.57911
5 5 2 -3526.02970
5 5 3 17340.16294
5 5 4 -53246.66952
5 5 5 106338.89246
5 5 6 -139773.52823
5 5 7 120269.52289
5 5 8 -67191.30684
5 5 9 24988.18495
5 5 10 -6761.78806
5 5 11 1175.12760")
Bmat <- data.matrix(read.table(text=Bijk, header=TRUE))
b <- Bmat[ , "B"]
i <- Bmat[ , "i"]
j <- Bmat[ , "j"]
k <- Bmat[ , "k"]
sum(b * r1^(i-j) * ctr[1]^j * sqrt(lambda[2])^k * sqrt(lambda[1])^(-j-k+1), na.rm=TRUE)
} else {
warning("Zhang CEP is only available for 2D-data")
rep(NA_real_, length(CEPlevel))
}
if(any(CEPlevel != 0.5)) {
warning("Zhang CEP is only available for CEPlevel 0.5")
}
Zhang <- ifelse(CEPlevel == 0.5, Zhang50, NA_real_)
setNames(Zhang, CEPlevel)
}
#####-----------------------------------------------------------------------
## Krempasky CEP estimate from Krempasky, 2003
CEPKrempasky <-
function(CEPlevel=0.5, sigma=diag(2), accuracy=FALSE) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
## only available for 2D 50% POA=POI case
Krempasky50 <- if((p == 2L) && !accuracy) {
## estimated correlation, covariance and standard deviations
rho <- cov2cor(sigma)[1, 2]
covXY <- sigma[1, 2]
sigmaX <- sqrt(sigma[1, 1])
sigmaY <- sqrt(sigma[2, 2])
## rotation angle gamma and rotation matrix
gamma <- atan((-2*covXY + sqrt(4*covXY^2 + (sigmaY^2 - sigmaX^2)^2)) /
(sigmaY^2 - sigmaX^2))
A <- cbind(c(cos(gamma), sin(gamma)), c(-sin(gamma), cos(gamma)))
## covariance matrix, correlation and standard deviations of rotated data
sigmaRot <- t(A) %*% sigma %*% A # covariance matrix
rhoPrime <- cov2cor(sigmaRot)[1, 2] # correlation
sigmaPrime <- sqrt(sigmaRot[1, 1])
# isTRUE(all.equal(sigmaRot[1, 1], sigmaRot[2, 2])) # supposed to be equal
CEP00 <- sigmaPrime*1.1774100225154746910115693264596996377473856893858
C2 <- 0.3267132048600136726456919696354558579811249664099
C4 <- 0.0568534980324428522403361717417556935145611584613
## CEP00 <- sigmaPrime*sqrt(2*log(2)) # sigmaPrime * qRayleigh(0.5, 1)
## C2 <- 0.5*(1 - log(2)/2)
## C4 <- C2*(log(2) - log(2)^2/4 - 0.5) + 3/8 - (9/16)*log(2) + (3/16)*log(2)^2 - (1/64)*log(2)^3 - C2^2*log(2)/2
CEP00*(1 - 0.5*C2*rhoPrime^2 - 0.5*(C4 + 0.25*C2^2)*rhoPrime^4)
} else {
if(p != 2L) {
warning("Krempasky CEP is only available for 2D-data")
}
if(accuracy) {
warning("Krempasky CEP is only available for accuracy=FALSE")
}
rep(NA_real_, length(CEPlevel))
}
if(any(CEPlevel != 0.5)) {
warning("Krempasky CEP is only available for CEPlevel 0.5")
}
Krempasky <- ifelse(CEPlevel == 0.5, Krempasky50, NA_real_)
setNames(Krempasky, CEPlevel)
}
#####-----------------------------------------------------------------------
## Ignani estimate from Ignani (2010)
CEPIgnani <-
function(CEPlevel=0.5, sigma=diag(2), accuracy=FALSE) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
## make sure eigenvalues >= 0 when very small
ev <- eigen(sigma)$values # eigenvalues
lambda <- ev*sign(ev)
alpha <- sqrt(lambda[2] / lambda[1])
beta <- if(p == 3L) { sqrt(lambda[3] / lambda[2]) } else { 0 }
betaVec <- c(1, beta, beta^2, beta^3)
## coefficients for polynomials in table 1
tbl1 <- c("coef R0.5 R0.9 R0.95 R0.99
c11 0.6754 1.6494 1.9626 2.5686
c12 -0.1547 0.0332 -0.0906 -0.1150
c13 0.2616 1.3376 1.3214 -0.3475
c14 1.0489 -0.8445 -0.3994 1.3570
c21 -0.0208 -0.0588 0.0100 0.1479
c22 1.1739 -0.5605 0.2722 0.9950
c23 1.9540 -4.7170 -5.4821 1.3223
c24 -5.5678 5.7135 3.9732 -4.8917
c31 1.1009 0.3996 0.0700 -0.4285
c32 -2.6375 1.5739 0.0462 -1.9795
c33 -1.4838 5.3623 7.1658 -1.1104
c34 6.5837 -7.9347 -5.7194 6.9617
c41 -0.5821 0.1636 0.4092 0.7371
c42 1.5856 -1.0747 -0.1953 1.2495
c43 -0.0678 -1.7785 -3.0134 -0.2061
c44 -2.3324 3.2388 2.4661 -2.8968")
coefDF <- read.table(text=tbl1, header=TRUE)
getR <- function(level) {
column <- paste0("R", level)
c1 <- coefDF[ 1:4, column]
c2 <- coefDF[ 5:8, column]
c3 <- coefDF[ 9:12, column]
c4 <- coefDF[13:16, column]
sqrt(lambda[1]) * (crossprod(c1, betaVec) +
crossprod(c2, betaVec)*alpha +
crossprod(c3, betaVec)*alpha^2 +
crossprod(c4, betaVec)*alpha^3)
}
Ignani <- if((p %in% c(2L, 3L)) && !accuracy) {
vapply(CEPlevel, function(x) {
if(x %in% c(0.5, 0.9, 0.95, 0.99)) { getR(x) } else { NA_real_ }
}, numeric(1))
} else {
if(!(p %in% c(2L, 3L))) {
warning("Ignani CEP is only available for 2D/3D-data")
}
if(accuracy) {
warning("Ignani CEP is only available for accuracy=FALSE")
}
rep(NA_real_, length(CEPlevel))
}
if(!any(CEPlevel %in% c(0.5, 0.9, 0.95, 0.99))) {
warning("Ignani CEP is only available for CEPlevel 0.5, 0.9, 0.95, 0.99")
}
setNames(Ignani, CEPlevel)
}
#####-----------------------------------------------------------------------
## RMSE-based estimate from van Diggelen (2007)
## CEP = sqrt(qchisq(0.5, 2)) * RMSEx
## RMSExy = sqrt(2) * RMSEx = sqrt(2) * (1/sqrt(qchisq(0.5, 2)))*CEP
CEPRMSE <-
function(CEPlevel=0.5, sigma=diag(2), accuracy=FALSE, xy) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
## root mean squared error
RMSE_AT <- sqrt(sum(colMeans(xy^2))) # non-centered data
RMSE_AF <- sqrt(sum(diag(sigma))) # centered data
## conversion factor from RMSExy or RMSExyz to CEP
## 2D: (1/sqrt(2)) * qRayleigh(CEPlevel, scale=1)
## for 50%, this is just sqrt(log(2))
## 3D: (1/sqrt(3)) * qMaxwell(CEPlevel, sigma=1)
RMSE2CEP <- (1/sqrt(p)) * sqrt(qchisq(CEPlevel, df=p))
RMSE <- if(accuracy) {
## POA != POI -> essentially Rice (2D) / offset sphere (3D)
RMSE2CEP * RMSE_AT
} else {
## POA = POI -> essentially Rayleigh (2D) / Maxwell-Boltzmann(3D)
RMSE2CEP * RMSE_AF
}
setNames(RMSE, CEPlevel)
}
#####-----------------------------------------------------------------------
## Ethridge CEP estimate from Ethridge (1983) after Puhek (1992)
CEPEthridge <-
function(CEPlevel=0.5, accuracy=FALSE, xy) {
lnDTC <- if(accuracy) { # log distance to group center (radius)
rSqSum <- sqrt(rowSums(xy^2)) # log radii to origin = point of aim
log(rSqSum)
} else {
log(getDistToCtr(xy)) # log radii to group center
}
mLnDTC <- mean(lnDTC) # mean log radius
medLnDTC <- median(lnDTC) # median log radius
varLnDTC <- var(lnDTC) # variance log radius
## weighted mean after Hogg (1967)
## sample kurtosis log radius
kLnDTC <- mean((lnDTC - mLnDTC)^4) / mean((lnDTC - mLnDTC)^2)^2
dHogg <- pmax(1 + (0.03 * (kLnDTC-3)^3 * (lnDTC-medLnDTC)^2 / varLnDTC), 0.01)
wHogg <- (1/dHogg) / sum(1/dHogg) # weighting factors
uHogg <- sum(wHogg * lnDTC) # log median radius estimate
Ethridge50 <- exp(uHogg)
if(any(CEPlevel != 0.5)) {
warning("Ethridge CEP estimate is only available for CEPlevel 0.5")
}
Ethridge <- ifelse(CEPlevel == 0.5, Ethridge50, NA_real_)
setNames(Ethridge, CEPlevel)
}
#####-----------------------------------------------------------------------
## modified RAND-234 CEP estimate for 50% from Williams, 1997
## using the semi-major and semi-minor axes of the error ellipse (PCA)
CEPRAND <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
## make sure eigenvalues >= 0 when very small
ev <- eigen(sigma)$values # eigenvalues
lambda <- ev*sign(ev)
RAND50 <- if(p == 2L) { # only available for 2D case
RAND50MPI <- 0.5620*sqrt(lambda[1]) + 0.6152*sqrt(lambda[2])
if(accuracy) { # take systematic location bias into account
bias <- sqrt(sum(ctr^2)) / RAND50MPI
if(bias > 2.2) {
warning(c("RAND location bias estimate is ",
round(bias, 2), " (> 2.2),\n",
"more than what RAND CEP should be considered for"))
}
## cubic regression to take bias into account
RAND50MPI * (1.0039 - 0.0528*bias + 0.4786*bias^2 - 0.0793*bias^3)
} else { # ignore location bias
RAND50MPI
} # if(accuracy)
} else { # 1D/3D case
warning("RAND CEP estimate is only available for 2D-data")
rep(NA_real_, length(CEPlevel))
}
if(any(CEPlevel != 0.5)) {
warning("RAND CEP estimate is only available for CEPlevel 0.5")
}
RAND <- ifelse(CEPlevel == 0.5, RAND50, NA_real_)
setNames(RAND, CEPlevel)
}
#####-----------------------------------------------------------------------
## Valstar CEP estimate for 50% from Williams, 1997
## using the semi-major and semi-minor axes of the error ellipse (PCA)
CEPValstar <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
aspRat <- sqrt(kappa(sigma, exact=TRUE))
## make sure eigenvalues >= 0 when very small
ev <- eigen(sigma)$values # eigenvalues
lambda <- ev*sign(ev)
Valstar50 <- if(p == 2L) { # only available for 2D case
ValstarMPI <- if((1/aspRat) <= 0.369) {
0.675*sqrt(lambda[1]) + sqrt(lambda[2])/(1.2*sqrt(lambda[1]))
} else {
0.5620*sqrt(lambda[1]) + 0.6152*sqrt(lambda[2]) # almost RAND
}
if(accuracy) { # POA != POI
Valstar <- sqrt(ValstarMPI^2 + sum(ctr^2))
} else { # POA = POI
ValstarMPI
} # if(accuracy)
} else { # 1D/3D case
warning("Valstar CEP estimate is only available for 2D-data")
rep(NA_real_, length(CEPlevel))
}
if(any(CEPlevel != 0.5)) {
warning("Valstar CEP estimate is only available for CEPlevel 0.5")
}
Valstar <- ifelse(CEPlevel == 0.5, Valstar50, NA_real_)
setNames(Valstar, CEPlevel)
}
# ## Siouris, GM. 1993. Appendix A
# CEPSiouris <- function(sigma) {
# sigma <- as.matrix(sigma)
# p <- ncol(sigma)
# aspRat <- sqrt(kappa(sigma, exact=TRUE))
#
# ## make sure eigenvalues >= 0 when very small
# ev <- eigen(sigma)$values # eigenvalues
# lambda <- ev*sign(ev)
#
# if(p == 2L) { # 2D
# theta <- 0.5 * atan(2*sqrt(lambda[1]*lambda[2]) / (lambda[1] - lambda[2]))
# 0.589*(sqrt(lambda[1]*cos(theta)^2 + lambda[2]*sin(theta)^2) +
# sqrt(lambda[1]*sin(theta)^2 + lambda[2]*cos(theta)^2))
# } else if(p == 3L) { # 3D
# varTotal <- sum(lambda)
# V <- 2*sum(lambda^2) / varTotal^2
# sqrt(varTotal*(1-(V/9)^3))
# }
# }
#
# ## Shultz, ME. 1963. Circular error probability of a quantity affected by a bias
# CEPShultz <- function(ctr, sigma) {
# sigma <- as.matrix(sigma)
# p <- ncol(sigma)
# aspRat <- sqrt(kappa(sigma, exact=TRUE))
#
# ## make sure eigenvalues >= 0 when very small
# ev <- eigen(sigma)$values # eigenvalues
# lambda <- ev*sign(ev)
#
# muH <- sqrt(sum(ctr^2))
# sdC <- mean(sqrt(lambda))
# CE90 <- 2.1272*sdC + 0.1674*muH + 0.3623*(muH^2/sdC) - 0.055*(muH^3/sdC)
# }
#
# ## Ager, TP. 2004. An Analysis of Metric Accuracy Definitions
# ## and Methods of Computation NIMA InnoVision
# CEPAger <- function(ctr, sigma) {
# sigma <- as.matrix(sigma)
# p <- ncol(sigma)
# aspRat <- sqrt(kappa(sigma, exact=TRUE))
#
# ## make sure eigenvalues >= 0 when very small
# ev <- eigen(sigma)$values # eigenvalues
# lambda <- ev*sign(ev)
#
# muH <- sqrt(sum(ctr^2))
# sdC <- mean(sqrt(lambda))
#
# CE90 <- if(muH/sdC <= 0.1) {
# 2.1460*sdC
# } else if(muH/sdC <= 3) {
# 2.1272*sdC + 0.1674*muH + 0.3623*(muH^2/sdC) - 0.055*(muH^3/sdC)
# } else {
# 0.9860*muH + 1.4548*sdC
# }
# }
dwoll/shotGroups documentation built on Feb. 16, 2024, 2:21 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions CEPValstar CEPRAND CEPEthridge CEPRMSE CEPIgnani CEPKrempasky CEPZhang CEPRayleigh CEPGrubbsLiu CEPGrubbsPearson CEPGrubbsPatnaik CEPCorrNormal
#####-----------------------------------------------------------------------
## CEP estimate based on correlated bivariate normal distribution - mvnEll.R, hoyt.R
CEPCorrNormal <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
CorrNormal <- if(accuracy) { # POA != POI
## quantile from offset circle probability - most general case
qmvnEll(CEPlevel, mu=numeric(p), sigma=sigma, e=diag(p), x0=ctr)
} else { # POA == POI
if(p == 2L) { # 2D -> exact Hoyt distribution
HP <- getHoytParam(sigma)
qHoyt(CEPlevel, qpar=HP$q, omega=HP$omega)
} else { # 1D/3D case
qmvnEll(CEPlevel, mu=numeric(p), sigma=sigma, e=diag(p), x0=numeric(p))
}
}
setNames(CorrNormal, CEPlevel)
}
####-----------------------------------------------------------------------
## Grubbs estimates (Grubbs, 1964, p54, p55-56) - grubbs.R
## Grubbs-Patnaik CEP estimate based on Patnaik two-moment central chi^2 approximation
CEPGrubbsPatnaik <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE) {
GPP <- getGrubbsParam(sigma, ctr=ctr, accuracy=accuracy)
GrubbsPatnaik <- qChisqGrubbs(CEPlevel, m=GPP$m, v=GPP$v, n=GPP$n, type="Patnaik")
setNames(GrubbsPatnaik, CEPlevel)
}
## Grubbs-Pearson CEP estimate based on Pearson three-moment central chi^2 approximation
CEPGrubbsPearson <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE) {
GPP <- getGrubbsParam(sigma, ctr=ctr, accuracy=accuracy)
GrubbsPearson <- qChisqGrubbs(CEPlevel, m=GPP$m, v=GPP$v, nPrime=GPP$nPrime, type="Pearson")
setNames(GrubbsPearson, CEPlevel)
}
## Grubbs-Liu CEP estimate based on four-moment non-central chi^2
## approximation (Liu, Tang & Zhang, 2009)
CEPGrubbsLiu <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE) {
GPP <- getGrubbsParam(sigma, ctr=ctr, accuracy=accuracy)
GrubbsLiu <- qChisqGrubbs(CEPlevel, m=GPP$m, v=GPP$v, muX=GPP$muX,
varX=GPP$varX, l=GPP$l, delta=GPP$delta, type="Liu")
setNames(GrubbsLiu, CEPlevel)
}
#####-----------------------------------------------------------------------
## Rayleigh CEP estimate from Singh, 1992 - rayleigh.R, rice.R, maxwell.R
CEPRayleigh <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE, doRob=FALSE, xy) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
Rayleigh <- if(p == 1L) { # 1D
if(accuracy) { # POA != POI -> non-central F distribution
RiceParam <- getRiceParam(xy, doRob=doRob)
N <- length(xy)
NU <- RiceParam$nu
SD <- setNames(RiceParam$sigma["sigma"], NULL)
NCP <- 1*(NU/SD)^2 # N*f^2, N=1 because interested in individual shot, not mean
SD * sqrt(qf(CEPlevel, df1=1, df2=N-1, ncp=NCP))
} else { # POA = POI -> half normal / central F distribution
N <- length(xy)
SD <- getRayParam(xy, doRob=doRob)$sigma["sigma"]
SD * sqrt(qchisq(CEPlevel, df=p))
}
} else if(p == 2L) { # 2D
if(accuracy) { # POA != POI -> Rice distribution
RiceParam <- getRiceParam(xy, doRob=doRob)
qRice(CEPlevel, nu=RiceParam$nu, sigma=RiceParam$sigma["sigma"])
} else { # POA = POI -> Rayleigh
RayParam <- getRayParam(xy, doRob=doRob)
qRayleigh(CEPlevel, scale=RayParam$sigma["sigma"])
}
} else if(p == 3L) { # 3D
MaxParam <- getRayParam(xy, doRob=doRob)
if(accuracy) { # POA != POI -> offset sphere probability
## circular covariance matrix with estimated M-B param sigma
sigMat <- diag(rep(MaxParam$sigma["sigma"]^2, p))
qmvnEll(CEPlevel, sigma=sigMat, mu=numeric(p), x0=ctr, e=diag(p))
} else { # POA = POI -> Maxwell-Boltzmann distribution
qMaxwell(CEPlevel, sigma=MaxParam$sigma["sigma"])
}
} else {
warning("Rayleigh CEP is only available for 1D/2D/3D-data")
rep(NA_real_, length(CEPlevel))
}
setNames(Rayleigh, CEPlevel)
}
#####-----------------------------------------------------------------------
## Zhang estimate from Zhang & An (2012) p570 eq. 7
CEPZhang <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE,
isotropic=FALSE, doRob=FALSE, xy) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
## make sure eigenvalues >= 0 when very small
ev <- eigen(sigma)$values # eigenvalues
lambda <- ev*sign(ev)
if(!accuracy) {
ctr <- c(0, 0)
}
Zhang50 <- if((p == 2L) && isotropic) { # assume equal variances
sigma <- getRayParam(xy=xy, level=CEPlevel, doRob=doRob)[["sigma"]]["sigma"]
r0 <- sqrt(sum((ctr/sigma)^2)) # p 570
print(r0)
if(r0 <= sqrt(5)) {
A0 <- 1.17712063
A1 <- 0.01635468
A2 <- 0.21387229
A3 <- 0.13972335
A4 <- -0.08426949
A5 <- 0.01284512
sigma*(A0 + A1*r0 + A2*r0^2 + A3*r0^3 + A4*r0^4 + A5*r0^5)
} else {
if(r0 <= 3) {
d <- 4
B0 <- 2.2458024
B1 <- 0.2200099
B2 <- -0.01003873
B3 <- 0.00078179
B4 <- -0.00005397
B5 <- 0.00000204
} else if(r0 <= 4) {
d <- 9
B0 <- 3.165246
B1 <- 0.15761986
B2 <- -0.00387106
B3 <- 0.00017928
B4 <- -0.00000822
B5 <- 0.00000022
} else if(r0 <= 5) {
d <- 16
B0 <- 4.1243779
B1 <- 0.12114747
B2 <- -0.00177091
B3 <- 0.00005016
B4 <- -0.0000015
B5 <- 0.00000003
} else if(r0 <= sqrt(45)) {
d <- 25
B0 <- 5.099676
B1 <- 0.09800860
B2 <- -0.00093530
B3 <- 0.00001679
B4 <- -0.00000029
B5 <- 0
} else {
warning("Isotropic case only defined up to r0 <= sqrt(45)")
d <- NA_real_
B0 <- NA_real_
B1 <- NA_real_
B2 <- NA_real_
B3 <- NA_real_
B4 <- NA_real_
B5 <- NA_real_
}
sigma * (B0 + B1*(r0^2 - d) +
B2*(r0^2 - d)^2 +
B3*(r0^2 - d)^3 +
B4*(r0^2 - d)^4 +
B5*(r0^2 - d)^5)
}
} else if((p == 2L)) { # do not assume equal variances
r1 <- sqrt(sum(ctr^2/lambda)) # p568
## coefficients for eq. 7 (personal communication)
Bijk <- c("i j k B
0 0 0 -1.50800
0 0 1 37.06914
0 0 2 -261.85730
0 0 3 992.70823
0 0 4 -2127.35633
0 0 5 2432.20862
0 0 6 -866.72077
0 0 7 -1043.20376
0 0 8 1007.57426
0 0 9 201.09441
0 0 10 -547.24802
0 0 11 178.38609
0 1 0 30.73407
0 1 1 -559.72454
0 1 2 4340.04504
0 1 3 -18573.68822
0 1 4 48237.39491
0 1 5 -78643.98158
0 1 6 80365.58824
0 1 7 -51299.58759
0 1 8 24063.13327
0 1 9 -14082.77353
0 1 10 8340.22172
0 1 11 -2217.07894
0 2 0 -140.53635
0 2 1 2615.24869
0 2 2 -21136.32531
0 2 3 95616.46894
0 2 4 -266698.90970
0 2 5 478107.64088
0 2 6 -557690.04739
0 2 7 425021.90498
0 2 8 -221989.27415
0 2 9 95629.50576
0 2 10 -38007.28570
0 2 11 8670.58029
0 3 0 307.16842
0 3 1 -5592.48090
0 3 2 45454.19525
0 3 3 -210124.94755
0 3 4 605819.65479
0 3 5 -1134958.18022
0 3 6 1399681.23997
0 3 7 -1135207.32583
0 3 8 611593.91350
0 3 9 -236561.30418
0 3 10 74259.96471
0 3 11 -14670.14344
0 4 0 -338.98709
0 4 1 5886.07371
0 4 2 -46868.89957
0 4 3 216011.08369
0 4 4 -628036.11139
0 4 5 1196820.36257
0 4 6 -1511973.91364
0 4 7 1259646.73040
0 4 8 -686657.86707
0 4 9 251291.10861
0 4 10 -67419.94690
0 4 11 11638.94476
0 5 0 150.64942
0 5 1 -2481.96624
0 5 2 19094.1026
0 5 3 -86299.71627
0 5 4 248634.95180
0 5 5 -473172.89850
0 5 6 600573.05422
0 5 7 -504295.29927
0 5 8 275477.36039
0 5 9 -97818.1152801
0 5 10 23834.43367
0 5 11 -3696.11403
1 0 0 7.30917
1 0 1 -121.95931
1 0 2 840.42710
1 0 3 -3062.11230
1 0 4 6186.79357
1 0 5 -6054.59489
1 0 6 -208.97884
1 0 7 6535.71511
1 0 8 -5133.24462
1 0 9 -260.95477
1 0 10 1917.57177
1 0 11 -645.78090
1 1 0 -110.31830
1 1 1 2017.70282
1 1 2 -15683.19635
1 1 3 67488.60741
1 1 4 -176701.28435
1 1 5 290580.88672
1 1 6 -298955.31158
1 1 7 191453.37431
1 1 8 -90293.34635
1 1 9 54007.00700
1 1 10 -32486.81832
1 1 11 8681.02815
1 2 0 498.12492
1 2 1 -9427.49801
1 2 2 77155.89222
1 2 3 -354370.79904
1 2 4 10006477.27304
1 2 5 -1840218.83052
1 2 6 2190104.05741
1 2 7 -1701102.53064
1 2 8 899284.68715
1 2 9 -383894.11296
1 2 10 149222.77548
1 2 11 -33723.07562
1 3 0 -1035.26688
1 3 1 19430.28736
1 3 2 -161868.38075
1 3 3 767936.17634
1 3 4 -2276698.65746
1 3 5 4389361.94030
1 3 6 -5567199.26107
1 3 7 4629872.46270
1 3 8 -2529969.48968
1 3 9 961497.85815
1 3 10 -285581.48386
1 3 11 54243.94576
1 4 0 1101.18523
1 4 1 -19807.40970
1 4 2 162886.78484
1 4 3 -776205.30786
1 4 4 2336828.66380
1 4 5 -4612153.15252
1 4 6 6025709.64684
1 4 7 -5168691.92531
1 4 8 2862608.29402
1 4 9 -1024135.74071
1 4 10 251524.75518
1 4 11 -39658.01666
1 5 0 -487.75337
1 5 1 8290.53777
1 5 2 -65870.58342
1 5 3 308197.20491
1 5 4 -921057.63790
1 5 5 1819395.02951
1 5 6 -2393829.33582
1 5 7 2073889.98097
1 5 8 -1152239.51773
1 5 9 398819.42046
1 5 10 -86623.36549
1 5 11 11513.64407
2 0 0 -8.39680
2 0 1 137.59442
2 0 2 -926.41844
2 0 3 3265.94858
2 0 4 -6100.14231
2 0 5 4278.05349
2 0 6 4933.42238
2 0 7 -12641.00423
2 0 8 8961.99007
2 0 9 -256.37747
2 0 10 -2517.52899
2 0 11 872.69812
2 1 0 151.85852
2 1 1 -2764.43629
2 1 2 21418.68153
2 1 3 -92079.73132
2 1 4 240640.93680
2 1 5 -392722.30429
2 1 6 395837.40643
2 1 7 -243259.82654
2 1 8 112039.90079
2 1 9 -74462.27711
2 1 10 48378.10091
2 1 11 -13174.53964
2 2 0 -696.05269
2 2 1 13235.79923
2 2 2 -108719.94280
2 2 3 502418.93939
2 2 4 -1437223.47953
2 2 5 2643175.11851
2 2 6 -3154239.38546
2 2 7 2448763.15222
2 2 8 -1297826.56619
2 2 9 566806.86704
2 2 10 -228178.84490
2 2 11 52471.09248
2 3 0 1400.85858
2 3 1 -26862.13543
2 3 2 227395.73358
2 3 3 -1096233.13091
2 3 4 3302421.47593
2 3 5 -6460706.72431
2 3 6 8294685.64827
2 3 7 -6962294.81091
2 3 8 3829716.07859
2 3 9 -1463576.88726
2 3 10 437814.27491
2 3 11 -83738.96835
2 4 0 -1432.43625
2 4 1 26735.89822
2 4 2 -226271.58166
2 4 3 1106145.22484
2 4 4 -3409574.42966
2 4 5 6873381.61298
2 4 6 -9144047.77589
2 4 7 7956459.96580
2 4 8 -4444585.05705
2 4 9 1586114.68543
2 4 10 -381893.75700
2 4 11 58951.23793
2 5 0 622.91812
2 5 1 -11035.33184
2 5 2 90838.15831
2 5 3 -438500.09965
2 5 4 1348016.12322
2 5 5 -2731098.77184
2 5 6 3673259.51515
2 5 7 -3238879.5041
2 5 8 1817714.84110
2 5 9 -624450.71184
2 5 10 129495.99917
2 5 11 -15978.20880
3 0 0 4.80967
3 0 1 -77.99167
3 0 2 516.96171
3 0 3 -1760.05286
3 0 4 2978.06661
3 0 5 -938.44505
3 0 6 -5613.89495
3 0 7 10344.99174
3 0 8 -6992.49989
3 0 9 -507.98776
3 0 10 1593.73407
3 0 11 -563.14246
3 1 0 -106.25020
3 1 1 1916.56093
3 1 2 -14700.28712
3 1 3 762447.47246
3 1 4 -160293.41157
3 1 5 253462.21660
3 1 6 -240535.58640
3 1 7 131669.95187
3 1 8 -55767.81377
3 1 9 47198.42290
3 1 10 -35116.21231
3 1 11 9820.81246
3 2 0 511.67374
3 2 1 -9673.11429
3 2 2 78791.45853
3 2 3 -360766.93776
3 2 4 1020321.00332
3 2 5 -1846866.00449
3 2 6 2154711.78200
3 2 7 -1626738.07815
3 2 8 852404.97838
3 2 9 -395825.66688
3 2 10 175015.53634
3 2 11 -41872.09886
3 3 0 -1027.83063
3 3 1 19916.0826
3 3 2 -168649.08473
3 3 3 810228.93041
3 3 4 -2426699.32685
3 3 5 4707457.17127
3 3 6 -5975393.05967
3 3 7 4953340.57762
3 3 8 -2714680.08066
3 3 9 1072562.25276
3 3 10 -347415.58697
3 3 11 70336.92513
3 4 0 1011.50846
3 4 1 -19520.35600
3 4 2 167740.84261
3 4 3 -825219.08597
3 4 4 2547613.20509
3 4 5 -5126865.40135
3 4 6 6791764.92505
3 4 7 -5879408.35900
3 4 8 8282864.52405
3 4 9 -1196814.35705
3 4 10 308131.74625
3 4 11 -51281.82176
3 5 0 -419.17227
3 5 1 7799.23359
3 5 2 -66151.92153
3 5 3 325082.85879
3 5 4 -1009967.60753
3 5 5 2058208.63394
3 5 6 -2775472.43675
3 5 7 2449243.21250
3 5 8 -1377607.57999
3 5 9 479600.83760
3 5 10 -104225.66822
3 5 11 13904.43044
4 0 0 -1.41012
4 0 1 22.70521
4 0 2 -148.46582
4 0 3 489.74223
4 0 4 -747.89192
4 0 5 -89.38056
4 0 6 2297.32756
4 0 7 -3725.25391
4 0 8 2452.57493
4 0 9 -251.29178
4 0 10 -464.39654
4 0 11 165.46501
4 1 0 36.91482
4 1 1 -657.65046
4 1 2 4965.58653
4 1 3 -20654.80936
4 1 4 51363.19899
4 1 5 -76954.43695
4 1 6 65622.22571
4 1 7 -27839.71942
4 1 8 9047.66734
4 1 9 -13506.39268
4 1 10 12032.45580
4 1 11 -3452.84254
4 2 0 -189.23560
4 2 1 3534.90907
4 2 2 -28340.16312
4 2 3 127301.53392
4 2 4 -351465.57188
4 2 5 616163.17964
4 2 6 -688181.61871
4 2 7 492046.07078
4 2 8 -251414.70516
4 2 9 129307.34993
4 2 10 -65040.55204
4 2 11 16270.97692
4 3 0 389.49639
4 3 1 -7517.54072
4 3 2 62883.56972
4 3 3 -297233.52178
4 3 4 872933.18776
4 3 5 -1653991.63581
4 3 6 2041991.02768
4 3 7 -1644559.84876
4 3 8 891362.10361
4 3 9 -372834.21958
4 3 10 136205.14245
4 3 11 -29615.35099
4 4 0 -375.58749
4 4 1 7349.92108
4 4 2 -63036.60176
4 4 3 307356.35957
4 4 4 -936607.41424
4 4 5 1854763.70135
4 4 6 -2412369.0135
4 4 7 2051411.99056
4 4 8 -1138662.20737
4 4 9 431751.74719
4 4 10 -124649.03235
4 4 11 23056.76527
4 5 0 146.79700
4 5 1 -2827.54869
4 5 2 24329.22818
4 5 3 -120032.35962
4 5 4 372079.75852
4 5 5 -753454.52592
4 5 6 1006974.01091
4 5 7 -880553.11272
4 5 8 494280.92111
4 5 9 -176872.67948
4 5 10 42410.61903
4 5 11 -6478.29063
5 0 0 0.16415
5 0 1 -2.62846
5 0 2 17.01287
5 0 3 -54.90894
5 0 4 78.17490
5 0 5 32.16035
5 0 6 -297.35896
5 0 7 454.60877
5 0 8 -293.38420
5 0 9 -37.40186
5 0 10 44.46108
5 0 11 -15.65335
5 1 0 -4.95117
5 1 1 86.91975
5 1 2 -643.82509
5 1 3 2608.72592
5 1 4 -6236.22491
5 1 5 8724.68903
5 1 6 -6379.03883
5 1 7 1475.79282
5 1 8 45.93268
5 1 9 1381.76407
5 1 10 -1495.997
5 1 11 435.75596
5 2 0 27.39345
5 2 1 -503.58840
5 2 2 3956.27375
5 2 3 -17332.34425
5 2 4 46350.38848
5 2 5 -77842.78318
5 2 6 81803.47877
5 2 7 -53864.40054
5 2 8 26400.02294
5 2 9 -15729.73024
5 2 10 9100.16729
5 2 11 -2363.21625
5 3 0 -58.88322
5 3 1 1121.35889
5 3 2 -9198.10897
5 3 3 42452.68426
5 3 4 -121199.98901
5 3 5 221998.94268
5 3 6 -263183.30573
5 3 7 202902.59108
5 3 8 -108017.59066
5 3 9 48909.67004
5 3 10 -20479.54143
5 3 11 4749.49844
5 4 0 56.60322
5 4 1 -1105.18410
5 4 2 961.78170
5 4 3 -44835.28439
5 4 4 133645.11524
5 4 5 -257933.62666
5 4 6 325955.13839
5 4 7 -269490.82354
5 4 8 148047.67616
5 4 9 -59284.13086
5 4 10 19575.20456
5 4 11 -3990.41775
5 5 0 -20.77692
5 5 1 407.57911
5 5 2 -3526.02970
5 5 3 17340.16294
5 5 4 -53246.66952
5 5 5 106338.89246
5 5 6 -139773.52823
5 5 7 120269.52289
5 5 8 -67191.30684
5 5 9 24988.18495
5 5 10 -6761.78806
5 5 11 1175.12760")
Bmat <- data.matrix(read.table(text=Bijk, header=TRUE))
b <- Bmat[ , "B"]
i <- Bmat[ , "i"]
j <- Bmat[ , "j"]
k <- Bmat[ , "k"]
sum(b * r1^(i-j) * ctr[1]^j * sqrt(lambda[2])^k * sqrt(lambda[1])^(-j-k+1), na.rm=TRUE)
} else {
warning("Zhang CEP is only available for 2D-data")
rep(NA_real_, length(CEPlevel))
}
if(any(CEPlevel != 0.5)) {
warning("Zhang CEP is only available for CEPlevel 0.5")
}
Zhang <- ifelse(CEPlevel == 0.5, Zhang50, NA_real_)
setNames(Zhang, CEPlevel)
}
#####-----------------------------------------------------------------------
## Krempasky CEP estimate from Krempasky, 2003
CEPKrempasky <-
function(CEPlevel=0.5, sigma=diag(2), accuracy=FALSE) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
## only available for 2D 50% POA=POI case
Krempasky50 <- if((p == 2L) && !accuracy) {
## estimated correlation, covariance and standard deviations
rho <- cov2cor(sigma)[1, 2]
covXY <- sigma[1, 2]
sigmaX <- sqrt(sigma[1, 1])
sigmaY <- sqrt(sigma[2, 2])
## rotation angle gamma and rotation matrix
gamma <- atan((-2*covXY + sqrt(4*covXY^2 + (sigmaY^2 - sigmaX^2)^2)) /
(sigmaY^2 - sigmaX^2))
A <- cbind(c(cos(gamma), sin(gamma)), c(-sin(gamma), cos(gamma)))
## covariance matrix, correlation and standard deviations of rotated data
sigmaRot <- t(A) %*% sigma %*% A # covariance matrix
rhoPrime <- cov2cor(sigmaRot)[1, 2] # correlation
sigmaPrime <- sqrt(sigmaRot[1, 1])
# isTRUE(all.equal(sigmaRot[1, 1], sigmaRot[2, 2])) # supposed to be equal
CEP00 <- sigmaPrime*1.1774100225154746910115693264596996377473856893858
C2 <- 0.3267132048600136726456919696354558579811249664099
C4 <- 0.0568534980324428522403361717417556935145611584613
## CEP00 <- sigmaPrime*sqrt(2*log(2)) # sigmaPrime * qRayleigh(0.5, 1)
## C2 <- 0.5*(1 - log(2)/2)
## C4 <- C2*(log(2) - log(2)^2/4 - 0.5) + 3/8 - (9/16)*log(2) + (3/16)*log(2)^2 - (1/64)*log(2)^3 - C2^2*log(2)/2
CEP00*(1 - 0.5*C2*rhoPrime^2 - 0.5*(C4 + 0.25*C2^2)*rhoPrime^4)
} else {
if(p != 2L) {
warning("Krempasky CEP is only available for 2D-data")
}
if(accuracy) {
warning("Krempasky CEP is only available for accuracy=FALSE")
}
rep(NA_real_, length(CEPlevel))
}
if(any(CEPlevel != 0.5)) {
warning("Krempasky CEP is only available for CEPlevel 0.5")
}
Krempasky <- ifelse(CEPlevel == 0.5, Krempasky50, NA_real_)
setNames(Krempasky, CEPlevel)
}
#####-----------------------------------------------------------------------
## Ignani estimate from Ignani (2010)
CEPIgnani <-
function(CEPlevel=0.5, sigma=diag(2), accuracy=FALSE) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
## make sure eigenvalues >= 0 when very small
ev <- eigen(sigma)$values # eigenvalues
lambda <- ev*sign(ev)
alpha <- sqrt(lambda[2] / lambda[1])
beta <- if(p == 3L) { sqrt(lambda[3] / lambda[2]) } else { 0 }
betaVec <- c(1, beta, beta^2, beta^3)
## coefficients for polynomials in table 1
tbl1 <- c("coef R0.5 R0.9 R0.95 R0.99
c11 0.6754 1.6494 1.9626 2.5686
c12 -0.1547 0.0332 -0.0906 -0.1150
c13 0.2616 1.3376 1.3214 -0.3475
c14 1.0489 -0.8445 -0.3994 1.3570
c21 -0.0208 -0.0588 0.0100 0.1479
c22 1.1739 -0.5605 0.2722 0.9950
c23 1.9540 -4.7170 -5.4821 1.3223
c24 -5.5678 5.7135 3.9732 -4.8917
c31 1.1009 0.3996 0.0700 -0.4285
c32 -2.6375 1.5739 0.0462 -1.9795
c33 -1.4838 5.3623 7.1658 -1.1104
c34 6.5837 -7.9347 -5.7194 6.9617
c41 -0.5821 0.1636 0.4092 0.7371
c42 1.5856 -1.0747 -0.1953 1.2495
c43 -0.0678 -1.7785 -3.0134 -0.2061
c44 -2.3324 3.2388 2.4661 -2.8968")
coefDF <- read.table(text=tbl1, header=TRUE)
getR <- function(level) {
column <- paste0("R", level)
c1 <- coefDF[ 1:4, column]
c2 <- coefDF[ 5:8, column]
c3 <- coefDF[ 9:12, column]
c4 <- coefDF[13:16, column]
sqrt(lambda[1]) * (crossprod(c1, betaVec) +
crossprod(c2, betaVec)*alpha +
crossprod(c3, betaVec)*alpha^2 +
crossprod(c4, betaVec)*alpha^3)
}
Ignani <- if((p %in% c(2L, 3L)) && !accuracy) {
vapply(CEPlevel, function(x) {
if(x %in% c(0.5, 0.9, 0.95, 0.99)) { getR(x) } else { NA_real_ }
}, numeric(1))
} else {
if(!(p %in% c(2L, 3L))) {
warning("Ignani CEP is only available for 2D/3D-data")
}
if(accuracy) {
warning("Ignani CEP is only available for accuracy=FALSE")
}
rep(NA_real_, length(CEPlevel))
}
if(!any(CEPlevel %in% c(0.5, 0.9, 0.95, 0.99))) {
warning("Ignani CEP is only available for CEPlevel 0.5, 0.9, 0.95, 0.99")
}
setNames(Ignani, CEPlevel)
}
#####-----------------------------------------------------------------------
## RMSE-based estimate from van Diggelen (2007)
## CEP = sqrt(qchisq(0.5, 2)) * RMSEx
## RMSExy = sqrt(2) * RMSEx = sqrt(2) * (1/sqrt(qchisq(0.5, 2)))*CEP
CEPRMSE <-
function(CEPlevel=0.5, sigma=diag(2), accuracy=FALSE, xy) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
## root mean squared error
RMSE_AT <- sqrt(sum(colMeans(xy^2))) # non-centered data
RMSE_AF <- sqrt(sum(diag(sigma))) # centered data
## conversion factor from RMSExy or RMSExyz to CEP
## 2D: (1/sqrt(2)) * qRayleigh(CEPlevel, scale=1)
## for 50%, this is just sqrt(log(2))
## 3D: (1/sqrt(3)) * qMaxwell(CEPlevel, sigma=1)
RMSE2CEP <- (1/sqrt(p)) * sqrt(qchisq(CEPlevel, df=p))
RMSE <- if(accuracy) {
## POA != POI -> essentially Rice (2D) / offset sphere (3D)
RMSE2CEP * RMSE_AT
} else {
## POA = POI -> essentially Rayleigh (2D) / Maxwell-Boltzmann(3D)
RMSE2CEP * RMSE_AF
}
setNames(RMSE, CEPlevel)
}
#####-----------------------------------------------------------------------
## Ethridge CEP estimate from Ethridge (1983) after Puhek (1992)
CEPEthridge <-
function(CEPlevel=0.5, accuracy=FALSE, xy) {
lnDTC <- if(accuracy) { # log distance to group center (radius)
rSqSum <- sqrt(rowSums(xy^2)) # log radii to origin = point of aim
log(rSqSum)
} else {
log(getDistToCtr(xy)) # log radii to group center
}
mLnDTC <- mean(lnDTC) # mean log radius
medLnDTC <- median(lnDTC) # median log radius
varLnDTC <- var(lnDTC) # variance log radius
## weighted mean after Hogg (1967)
## sample kurtosis log radius
kLnDTC <- mean((lnDTC - mLnDTC)^4) / mean((lnDTC - mLnDTC)^2)^2
dHogg <- pmax(1 + (0.03 * (kLnDTC-3)^3 * (lnDTC-medLnDTC)^2 / varLnDTC), 0.01)
wHogg <- (1/dHogg) / sum(1/dHogg) # weighting factors
uHogg <- sum(wHogg * lnDTC) # log median radius estimate
Ethridge50 <- exp(uHogg)
if(any(CEPlevel != 0.5)) {
warning("Ethridge CEP estimate is only available for CEPlevel 0.5")
}
Ethridge <- ifelse(CEPlevel == 0.5, Ethridge50, NA_real_)
setNames(Ethridge, CEPlevel)
}
#####-----------------------------------------------------------------------
## modified RAND-234 CEP estimate for 50% from Williams, 1997
## using the semi-major and semi-minor axes of the error ellipse (PCA)
CEPRAND <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
## make sure eigenvalues >= 0 when very small
ev <- eigen(sigma)$values # eigenvalues
lambda <- ev*sign(ev)
RAND50 <- if(p == 2L) { # only available for 2D case
RAND50MPI <- 0.5620*sqrt(lambda[1]) + 0.6152*sqrt(lambda[2])
if(accuracy) { # take systematic location bias into account
bias <- sqrt(sum(ctr^2)) / RAND50MPI
if(bias > 2.2) {
warning(c("RAND location bias estimate is ",
round(bias, 2), " (> 2.2),\n",
"more than what RAND CEP should be considered for"))
}
## cubic regression to take bias into account
RAND50MPI * (1.0039 - 0.0528*bias + 0.4786*bias^2 - 0.0793*bias^3)
} else { # ignore location bias
RAND50MPI
} # if(accuracy)
} else { # 1D/3D case
warning("RAND CEP estimate is only available for 2D-data")
rep(NA_real_, length(CEPlevel))
}
if(any(CEPlevel != 0.5)) {
warning("RAND CEP estimate is only available for CEPlevel 0.5")
}
RAND <- ifelse(CEPlevel == 0.5, RAND50, NA_real_)
setNames(RAND, CEPlevel)
}
#####-----------------------------------------------------------------------
## Valstar CEP estimate for 50% from Williams, 1997
## using the semi-major and semi-minor axes of the error ellipse (PCA)
CEPValstar <-
function(CEPlevel=0.5, ctr=c(0, 0), sigma=diag(length(ctr)), accuracy=FALSE) {
sigma <- as.matrix(sigma)
p <- ncol(sigma)
aspRat <- sqrt(kappa(sigma, exact=TRUE))
## make sure eigenvalues >= 0 when very small
ev <- eigen(sigma)$values # eigenvalues
lambda <- ev*sign(ev)
Valstar50 <- if(p == 2L) { # only available for 2D case
ValstarMPI <- if((1/aspRat) <= 0.369) {
0.675*sqrt(lambda[1]) + sqrt(lambda[2])/(1.2*sqrt(lambda[1]))
} else {
0.5620*sqrt(lambda[1]) + 0.6152*sqrt(lambda[2]) # almost RAND
}
if(accuracy) { # POA != POI
Valstar <- sqrt(ValstarMPI^2 + sum(ctr^2))
} else { # POA = POI
ValstarMPI
} # if(accuracy)
} else { # 1D/3D case
warning("Valstar CEP estimate is only available for 2D-data")
rep(NA_real_, length(CEPlevel))
}
if(any(CEPlevel != 0.5)) {
warning("Valstar CEP estimate is only available for CEPlevel 0.5")
}
Valstar <- ifelse(CEPlevel == 0.5, Valstar50, NA_real_)
setNames(Valstar, CEPlevel)
}
# ## Siouris, GM. 1993. Appendix A
# CEPSiouris <- function(sigma) {
# sigma <- as.matrix(sigma)
# p <- ncol(sigma)
# aspRat <- sqrt(kappa(sigma, exact=TRUE))
#
# ## make sure eigenvalues >= 0 when very small
# ev <- eigen(sigma)$values # eigenvalues
# lambda <- ev*sign(ev)
#
# if(p == 2L) { # 2D
# theta <- 0.5 * atan(2*sqrt(lambda[1]*lambda[2]) / (lambda[1] - lambda[2]))
# 0.589*(sqrt(lambda[1]*cos(theta)^2 + lambda[2]*sin(theta)^2) +
# sqrt(lambda[1]*sin(theta)^2 + lambda[2]*cos(theta)^2))
# } else if(p == 3L) { # 3D
# varTotal <- sum(lambda)
# V <- 2*sum(lambda^2) / varTotal^2
# sqrt(varTotal*(1-(V/9)^3))
# }
# }
#
# ## Shultz, ME. 1963. Circular error probability of a quantity affected by a bias
# CEPShultz <- function(ctr, sigma) {
# sigma <- as.matrix(sigma)
# p <- ncol(sigma)
# aspRat <- sqrt(kappa(sigma, exact=TRUE))
#
# ## make sure eigenvalues >= 0 when very small
# ev <- eigen(sigma)$values # eigenvalues
# lambda <- ev*sign(ev)
#
# muH <- sqrt(sum(ctr^2))
# sdC <- mean(sqrt(lambda))
# CE90 <- 2.1272*sdC + 0.1674*muH + 0.3623*(muH^2/sdC) - 0.055*(muH^3/sdC)
# }
#
# ## Ager, TP. 2004. An Analysis of Metric Accuracy Definitions
# ## and Methods of Computation NIMA InnoVision
# CEPAger <- function(ctr, sigma) {
# sigma <- as.matrix(sigma)
# p <- ncol(sigma)
# aspRat <- sqrt(kappa(sigma, exact=TRUE))
#
# ## make sure eigenvalues >= 0 when very small
# ev <- eigen(sigma)$values # eigenvalues
# lambda <- ev*sign(ev)
#
# muH <- sqrt(sum(ctr^2))
# sdC <- mean(sqrt(lambda))
#
# CE90 <- if(muH/sdC <= 0.1) {
# 2.1460*sdC
# } else if(muH/sdC <= 3) {
# 2.1272*sdC + 0.1674*muH + 0.3623*(muH^2/sdC) - 0.055*(muH^3/sdC)
# } else {
# 0.9860*muH + 1.4548*sdC
# }
# }
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.