r - 在 R 中使用 optim() 重现 Fisher 关于球形数据的书的结果

Question

我正在尝试从“球形数据的统计分析”中重现结果。
我想计算球面中位数(公式可以看http://www.jstor.org/stable/2345577，等式1，这里不知道怎么写才对)。
我使用本书的 B1 数据集:

lat1=c(-26.4,-32.2,-73.1,-80.2,-71.1,-58.7,-40.8,-14.9,-66.1,-1.8,-52.1,-77.3,-68.8,-68.4,
   -29.2,-78.5,-65.4,-49,-67,-56.7,-80.5,-77.7,-6.9,-59.4,-5.6,-62.6,-74.7,-65.3,-71.6,
   -23.3,-74.3,-81,-12.7,-75.4,-85.9,-84.8,-7.4,-29.8,-85.2,-53.1,-38.3,-72.7,-60.2,-63.4,
   -17.2,-81.6,-40.4,-53.6,-56.2,-75.1)

long1=c(324,163.7,51.9,140.5,267.2,32,28.1,266.3,144.3,256.2,83.2,182.1,110.4,142.2,246.3,222.6,247.7,
    65.6,282.6,56.2,108.4,266,19.1,281.7,107.4,105.3,120.2,286.6,106.4,96.5,90.2,170.9,199.4,118.6,
    63.7,74.9,93.8,72.8,113.2,51.5,146.8,103.1,33.2,154.8,89.9,295.6,41.0,59.1,35.6,70.7)

library('sphereplot')
B1=data.frame(long=long1,lat=lat1)
a=sph2car(B1$long,B1$lat)
x=a[,1]
y=a[,2]
z=a[,3]

我首先检查数据:

sqrt(x^2+y^2+z^2)

data1=data.frame(x,y,z)

median.direction <- function(par, data1) {
sum(acos(par[1]*data1[,1]+par[2]*data1[,2]+par[3]*data1[,3]))
}

median.direction2=optim(par=c(0,0,0), fn=median.direction, data1=data1)    
result1=car2sph(median.direction2$par[1],median.direction2$par[2],median.direction2$par[3])

result1

“对于例 5.1 的数据(集合 Bl)，球面中位数
方向是(纬度。78.9°，长。98.4°)。”

我不知道我的错误在哪里:

我必须在 sph2car 中使用 colatitude 吗？
optim 在警告时表现良好吗？

编辑 :

Answer 1

这里有几件事情正在发生。
首先，当数据集中的所有纬度都 < 0 时，很难看出中位纬度如何是 +79°。所以要么你的问题有错别字，要么教科书有错误。
其次，您的数据集中(或多或少)靠近其中一个极点。在这种情况下，您估计经度的能力本质上会受到损害。考虑所有数据都在纬度 -90° 的极端情况。那么中位纬度将恰好是 -90°，但我们对中位经度一无所知。所以你的优化问题是经度有一个“浅最小值”(在这个数据集中)。也就是说，有许多经度非常接近于最小化您的目标函数。这是一个问题，因为大多数优化器使用局部最小值——它们在接近初始估计的目标函数中寻找最小值。因此，您得到的答案将取决于您的起点。
第三，鉴于上述情况，您最好使用更强大的优化器 (IMO)。在下面的示例中，我使用 nloptr(...)来自 nloptr包裹。它使用起来有点困难，但给出的结果对初始估计不太敏感。
为了演示这个问题，下面的代码运行了 100 次最小化，每次都随机选择一个起点，并绘制数据和 100 个“最小值”。

library(sphereplot)
library(nloptr)
f <- function(par, data1) {
  sum(acos(par[1]*data1[,1]+par[2]*data1[,2]+par[3]*data1[,3]))
}
opts <- list(algorithm="NLOPT_GN_ISRES",xtol_rel=1.0e-6, maxeval=10000)
# set up the plot
rgl.sphgrid()
points3d(x,y,z, col="red",size=5)

set.seed(1)    # for reproducibility
# 100 initial estimates, randomly distributed on the sphere
N <- 100
xyz.init <- sph2car(long=sample(-180:180,N),lat=sample(-90:90,N))
get.median <- function(i) {
  md     <- nloptr(x0=xyz.init[i,],eval_f=f,
                   lb=c(-1,-1,-1), ub=c(1,1,1),
                   data1=data1, opts=opts)
  xyz    <- md$solution
  lines3d(c(0,xyz[1]),c(0,xyz[2]),c(0,xyz[3]),col="green",lwd=2)
  median <- car2sph(xyz[1],xyz[2],xyz[3])
  cat(".")     # cheap and dirty progress bar...
  return(median)
}  
sph.med  <- do.call(rbind,lapply(1:nrow(xyz.init),get.median))
colMeans(sph.med)
#       long        lat     radius 
#  92.314309 -77.361522   0.998315 

您可以看到优化为“中位数”创建了一个估计的信封(锥体)。所有这些估计的平均值非常接近书中的结果(纬度符号除外)。
值得注意的是，尽管使用最多 10,000 次迭代，但优化通常不会收敛!!