这是一个很好的问题。
我将遵循@amoeba 的建议depth::med()并使用with引导空间中位数method="Spatial"。但是，有一点复杂：med当有重复的数据点时不喜欢它，所以我们不能做一个简单的引导。相反，我将绘制一个引导样本，然后将每个点抖动一小部分 - 小于每个点中的最小距离x和y原始数据样本中的维度 - 在计算空间中位数之前。
最后，我将计算覆盖指定比例 (95%) 自举中位数和绘图的最小椭圆。
library(depth) # for med()
library(MASS) # for cov.rob()
library(cluster) # for ellipsoidhull()
# create data
set.seed(1)
df <- data.frame(x = rnorm(200, mean = 4, sd = 1.5),
y = rnorm(200, mean = 1.4, sd = 2.5))
# find minimum distances in each dimension for later jittering
foo <- outer(X=df$x,Y=df$x,FUN=function(xx,yy)abs(xx-yy))
delta.x <- min(foo[upper.tri(foo)])/2
foo <- outer(X=df$y,Y=df$y,FUN=function(xx,yy)abs(xx-yy))
delta.y <- min(foo[upper.tri(foo)])/2
# bootstrap spatial medians, using jittering
n.boot <- 1000
pb <- winProgressBar(max=n.boot)
boot.med <- matrix(NA,nrow=n.boot,ncol=2)
for ( ii in 1:n.boot ) {
setWinProgressBar(pb,ii,paste(ii,"of",n.boot))
index <- sample(1:nrow(df),nrow(df),replace=TRUE)
bar <- df[index,] +
data.frame(x=runif(nrow(df),-delta.x,delta.x),
y=runif(nrow(df),-delta.y,delta.y))
boot.med[ii,] <- med(bar,method="Spatial")$median
}
close(pb)
# specify confidence level
pp <- 0.95
# find smallest ellipse containing the specified proportion of bootstrapped medians
fit <- cov.rob(boot.med, quantile.used = ceiling(pp*n.boot), method = "mve")
best_ellipse <- ellipsoidhull( boot.med[fit$best,] )
plot(df)
points(boot.med,pch=19,col="grey",cex=0.5)
points(df)
lines(predict(best_ellipse), col="red")
legend("bottomright",bg="white",pch=c(21,19,NA),
col=c("black","grey","red"),pt.bg=c("white",NA,NA),lwd=c(0,0,1),
legend=c("Observations","Bootstrapped medians","Confidence ellipse"))
最后，注意二元空间中位数是渐近正态分布的（Brown, 1983, JRSS, Series B），所以我们也可以省去上面的“jittered bootstrap”，直接计算椭圆，相信n=200是“足够渐近的”。如果我在接下来的几天找到时间，我可能会编辑这篇文章以包含这个参数置信椭圆。
