机器算法验证 - 如何创建方位角 QQ 图来测试球面点数据集的旋转对称性？ - 吾爱随笔录

如何创建方位角 QQ 图来测试球面点数据集的旋转对称性？

机器算法验证矩阵循环统计

2022-04-13 18:00:48

我正在尝试做一些 QQ 图和关于平均方向的旋转对称性的柯伊伯检验。这些是单位向量、球形数据。我正在努力理解的是从我自己的方向旋转，其中 φ 是方位角方向（0 到 360），θ 是向上/向下倾斜方向（0 到 180），到旋转/修改的 φ' 和 θ' . 我将把这个帖子限制在 QQ 情节问题上。

这是一组 10 个方位角/φ 角的小示例：146.6 198.0 182.5 148.5 183.0 130.5 344.9 355.5 82.3 166.2

所以，首先是图解法：QQ图。图片，包括下面的#1，来自 Fisher、Lewis 和 Embleton 1987（参考下文）。我有 φ 值，但绘图需要旋转/修改的 φ' 值。

图形和形式测试

现在在下面的图像#2 中是进行旋转的方法，方程 3.8。我想我可以将 ψ ₀设置为零，所以我可以使用“更简单”的方程 3.9。

旋转方程

但是方程 3.8 或 3.9 需要方程 3.5 和 3.6，它们分别位于图像 #3 和 #4 的下方。我能够计算结果长度 (R) 和 x-、y- 和 z-hat，所以我有这些，这让我可以得到 φ-hat 和 φ-hat。所以，我不确定如何处理方程 3.9 或 3.9，或者如何旋转我的数字。有什么建议么？为图片的长篇道歉。

等式 3.5 等

方程 3.6

Fisher, NI, T. Lewis 和 BJJ Embleton。1987. 球形数据的统计分析。剑桥，剑桥大学出版社。

1个回答

您只想研究一组球面点相对于它们的球面平均值的方位角。最直接的解决方案是求解球面三角形，其中是北极。 $P_i$ $\bar P$ $(N,\bar P, P_i)$ $N$

和（与北方的角度）的共同纬度和。设为它们之间的角度：它只是相同两点的经度之差。方位角，正东为零并且逆时针方向的角度由下式确定 $P_i$ $\bar P$ $a$ $b$ $\gamma$

\arctan_{2} (\sin (b) \cos (a) - \cos (b) \sin (a) \cos (γ), \sin (a) \sin (γ))

$\arctan_2(\sin(b)\cos(a) - \cos(b)\sin(a)\cos(\gamma),\ \sin(a)\sin(\gamma))$

其中是平面中点的角度。（这应该是该公式的数值稳定版本，但我尚未对其进行广泛测试。） $\arctan_2(y,x)$ $(x,y)$

在此示例中，根据分布在南半球和西半球的（对称）Fisher-von Mises 分布生成了个点集中在南半球和东半球。结果分布不是对称的。 $100$ $50$

图显示了点和方位角的 QQ 图。

平均点显示为红色三角形。

到之间的方位角（以弧度表示）。和之间形成了更宽的方位角。QQ 图显然不均匀（否则它会靠近对角线虚线），反映了球面点分布的双峰性。 $0$ $1$ $3$ $5$

生成此示例的R代码可用于为任何数据生成方位 QQ 图。它假设球坐标以数组中的行的形式提供；相关行由“phi”和“theta”索引。

#
# Spherical triangle, two sides and included angle given.
# Returns the angle `alpha` opposite `a`, in radians between 0 and 2*pi.
#
SAS <- function(a, gamma, b) {
  atan2(sin(b)*cos(a) - cos(b)*sin(a)*cos(gamma), sin(a)*sin(gamma)) %% (2*pi)
}
#
# Cartesian coordinate conversion (for generating points).
#
xyz.to.spherical <- function(xyz) {
  xyz <- matrix(xyz, nrow=3)
  x <- xyz[1,]; y <- xyz[2,]; z <- xyz[3,]
  r2 <- x^2 + y^2
  rho <- sqrt(r2 + z^2)
  theta <- pi/2 - atan2(z, sqrt(r2))
  phi <- atan2(y, x)
  theta[x==0 && y==0] <- sign(z) * pi/2
  return (rbind(rho, theta, phi))
}
#
# Generate random points on the sphere.
#
library(MASS)
set.seed(17)
n.1 <- 100
n.2 <- 50
mu.1 <- c(0,-1,-1/4) * 2        # Center of first distribution
mu.2 <- c(1,1,-1/2) * 5         # Center of the second distribution
Sigma <- outer(1:3, 1:3, "==")  # Identity covariance matrix
xyz.1 <- t(mvrnorm(n.1, mu.1, Sigma)) # Each column is a point
xyz.2 <- t(mvrnorm(n.2, mu.2, Sigma))
xyz <- cbind(xyz.1, xyz.2)      # The Cartesian coordinates
rtf <- xyz.to.spherical(xyz)    # The spherical coordinates (also in columns)
#
# Compute the spherical mean and the azimuths relative to that mean.
#
mean.rtf <- xyz.to.spherical(rowMeans(xyz))
a <- SAS(rtf["theta",], rtf["phi",]-mean.rtf["phi",], mean.rtf["theta",])
#
# Plot the data and a QQ plot of the azimuths.
#
par(mfrow=c(1,2))
plot(c(-pi, pi), c(-1,1), type="n", 
     xlab="Phi", ylab="Cos(theta)", main="Sample points and their mean")
abline(h=0, col="Gray") # The Equator
abline(v=0, col="Gray") # The Prime Meridian
points(rtf["phi",], cos(rtf["theta", ]), col="#00000080")
points(mean.rtf["phi",], cos(mean.rtf["theta",]), bg="Red", pch=24, cex=1.25)

plot(c(0,1), c(0,2*pi), type="n",
     xlab="Quantile", ylab="Azimuth", main="Azimuthal QQ Plot")
abline(c(0, 2*pi), lty=3, lwd=2, col="Gray")
points(seq(0, 1, along.with=a), sort(a))

其它你可能感兴趣的问题