获取椭圆协方差图上的椭圆顶点（由 `car::ellipse` 创建）
Obtain vertices of the ellipse on an ellipse covariance plot (created by `car::ellipse`)
按照this post可以画出给定形状矩阵(A)的椭圆：
library(car)
A <- matrix(c(20.43, -8.59,-8.59, 24.03), nrow = 2)
ellipse(c(-0.05, 0.09), shape=A, radius=1.44, col="red", lty=2, asp = 1)
现在如何得到这个椭圆的major/minor（major/minor轴和椭圆的一对交点）顶点？
仍然非常不确定这是否真的能回答问题，但这是我的尝试：
首先，将椭圆的中心定义为向量，以备后用：
center<-c(x=-0.05, y=0.09)
绘制椭圆并得到 "points" 的矩阵，其值足够接近现实点：
tmp<-ellipse(c(-0.05, 0.09), shape=A, radius=1.44, segments=1e3, col="red", lty=2,add=FALSE)
用它创建一个 data.table 并计算每个点到中心的距离 (point_x - center_x)² + (point_y - center_y )²:
dt <- data.table(tmp)
dt[,dist:={dx=x-center[1];dy=y-center[2];dx*dx+dy*dy}]
按距离对顶点排序：
setorder(dt,dist)
获取最小和最大点数：
> tail(dt,2)
x y dist
1: 4.990415 -6.138039 64.29517
2: -5.110415 6.318039 64.29517
> head(dt,2)
x y dist
1: 4.045722 3.41267 27.89709
2: -4.165722 -3.23267 27.89709
不要添加太多段，否则前两个值将是非常接近而不是相反的两个点。
从视觉效果来看，这听起来不太准确：
出于实用目的，@Tensibai 的回答可能已经足够好了。只需为 segments 参数使用足够大的值，以便这些点可以很好地近似于真实顶点。
如果你想要更严格的东西，你可以求解沿椭圆的位置，maximises/minimises 到中心的距离，由角度参数化。由于形状矩阵的存在，这比仅采用 angle={0, pi/2, pi, 3pi/2} 更复杂。不过也不算太难:
# location along the ellipse
# linear algebra lifted from the code for ellipse()
ellipse.loc <- function(theta, center, shape, radius)
{
vert <- cbind(cos(theta), sin(theta))
Q <- chol(shape, pivot=TRUE)
ord <- order(attr(Q, "pivot"))
t(center + radius*t(vert %*% Q[, ord]))
}
# distance from this location on the ellipse to the center
ellipse.rad <- function(theta, center, shape, radius)
{
loc <- ellipse.loc(theta, center, shape, radius)
(loc[,1] - center[1])^2 + (loc[,2] - center[2])^2
}
# ellipse parameters
center <- c(-0.05, 0.09)
A <- matrix(c(20.43, -8.59, -8.59, 24.03), nrow=2)
radius <- 1.44
# solve for the maximum distance in one hemisphere (hemi-ellipse?)
t1 <- optimize(ellipse.rad, c(0, pi - 1e-5), center=center, shape=A, radius=radius, maximum=TRUE)$m
l1 <- ellipse.loc(t1, center, A, radius)
# solve for the minimum distance
t2 <- optimize(ellipse.rad, c(0, pi - 1e-5), center=center, shape=A, radius=radius)$m
l2 <- ellipse.loc(t2, center, A, radius)
# other points obtained by symmetry
t3 <- pi + t1
l3 <- ellipse.loc(t3, center, A, radius)
t4 <- pi + t2
l4 <- ellipse.loc(t4, center, A, radius)
# plot everything
MASS::eqscplot(center[1], center[2], xlim=c(-7, 7), ylim=c(-7, 7), xlab="", ylab="")
ellipse(center, A, radius, col="red", lty=2)
points(rbind(l1, l2, l3, l4), cex=2, col="blue", lwd=2)
我知道这个问题已经解决了，但实际上有一个超级优雅的解决方案，只需几行如下。这种计算是精确的，没有任何形式的数值优化。
## target covariance matrix
A <- matrix(c(20.43, -8.59,-8.59, 24.03), nrow = 2)
E <- eigen(A, symmetric = TRUE) ## symmetric eigen decomposition
U <- E[[2]] ## eigen vectors, i.e., rotation matrix
D <- sqrt(E[[1]]) ## root eigen values, i.e., scaling factor
r <- 1.44 ## radius of original circle
Z <- rbind(c(r, 0), c(0, r), c(-r, 0), c(0, -r)) ## original vertices on major / minor axes
Z <- tcrossprod(Z * rep(D, each = 4), U) ## transformed vertices on major / minor axes
# [,1] [,2]
#[1,] -5.055136 6.224212
#[2,] -4.099908 -3.329834
#[3,] 5.055136 -6.224212
#[4,] 4.099908 3.329834
C0 <- c(-0.05, 0.09) ## new centre
Z <- Z + rep(C0, each = 4) ## shift to new centre
# [,1] [,2]
#[1,] -5.105136 6.314212
#[2,] -4.149908 -3.239834
#[3,] 5.005136 -6.134212
#[4,] 4.049908 3.419834
为了解释背后的数学原理，我将分3步：
- 这个椭圆是从哪里来的？
- Cholesky 分解方法及其缺点。
- 特征分解方法及其自然解释。
这个椭圆是哪里来的？
在实践中，这个椭圆可以通过对单位圆的一些线性变换得到x ^ 2 + y ^ 2 = 1。
Cholesky分解法及其缺点
## initial circle
r <- 1.44
theta <- seq(0, 2 * pi, by = 0.01 * pi)
X <- r * cbind(cos(theta), sin(theta))
## target covariance matrix
A <- matrix(c(20.43, -8.59,-8.59, 24.03), nrow = 2)
R <- chol(A) ## Cholesky decomposition
X1 <- X %*% R ## linear transformation
Z <- rbind(c(r, 0), c(0, r), c(-r, 0), c(0, -r)) ## original vertices on major / minor axes
Z1 <- Z %*% R ## transformed coordinates
## different colour per quadrant
g <- floor(4 * (1:nrow(X) - 1) / nrow(X)) + 1
## draw ellipse
plot(X1, asp = 1, col = g)
points(Z1, cex = 1.5, pch = 21, bg = 5)
## draw circle
points(X, col = g, cex = 0.25)
points(Z, cex = 1.5, pch = 21, bg = 5)
## draw axes
abline(h = 0, lty = 3, col = "gray", lwd = 1.5)
abline(v = 0, lty = 3, col = "gray", lwd = 1.5)
我们看到线性变换矩阵 R 似乎没有自然解释。圆的原始顶点不映射到椭圆的顶点。
特征分解方法及其自然解释
## initial circle
r <- 1.44
theta <- seq(0, 2 * pi, by = 0.01 * pi)
X <- r * cbind(cos(theta), sin(theta))
## target covariance matrix
A <- matrix(c(20.43, -8.59,-8.59, 24.03), nrow = 2)
E <- eigen(A, symmetric = TRUE) ## symmetric eigen decomposition
U <- E[[2]] ## eigen vectors, i.e., rotation matrix
D <- sqrt(E[[1]]) ## root eigen values, i.e., scaling factor
r <- 1.44 ## radius of original circle
Z <- rbind(c(r, 0), c(0, r), c(-r, 0), c(0, -r)) ## original vertices on major / minor axes
## step 1: re-scaling
X1 <- X * rep(D, each = nrow(X)) ## anisotropic expansion to get an axes-aligned ellipse
Z1 <- Z * rep(D, each = 4L) ## vertices on axes
## step 2: rotation
Z2 <- tcrossprod(Z1, U) ## rotated vertices on major / minor axes
X2 <- tcrossprod(X1, U) ## rotated ellipse
## different colour per quadrant
g <- floor(4 * (1:nrow(X) - 1) / nrow(X)) + 1
## draw rotated ellipse and vertices
plot(X2, asp = 1, col = g)
points(Z2, cex = 1.5, pch = 21, bg = 5)
## draw axes-aligned ellipse and vertices
points(X1, col = g)
points(Z1, cex = 1.5, pch = 21, bg = 5)
## draw original circle
points(X, col = g, cex = 0.25)
points(Z, cex = 1.5, pch = 21, bg = 5)
## draw axes
abline(h = 0, lty = 3, col = "gray", lwd = 1.5)
abline(v = 0, lty = 3, col = "gray", lwd = 1.5)
## draw major / minor axes
segments(Z2[1,1], Z2[1,2], Z2[3,1], Z2[3,2], lty = 2, col = "gray", lwd = 1.5)
segments(Z2[2,1], Z2[2,2], Z2[4,1], Z2[4,2], lty = 2, col = "gray", lwd = 1.5)
在这里我们看到，在变换的两个阶段，顶点仍然映射到顶点。正是基于这样 属性 我们在一开始就给出了巧妙的解决方案。
按照this post可以画出给定形状矩阵(A)的椭圆：
library(car)
A <- matrix(c(20.43, -8.59,-8.59, 24.03), nrow = 2)
ellipse(c(-0.05, 0.09), shape=A, radius=1.44, col="red", lty=2, asp = 1)
现在如何得到这个椭圆的major/minor（major/minor轴和椭圆的一对交点）顶点？
仍然非常不确定这是否真的能回答问题，但这是我的尝试：
首先，将椭圆的中心定义为向量，以备后用：
center<-c(x=-0.05, y=0.09)
绘制椭圆并得到 "points" 的矩阵，其值足够接近现实点：
tmp<-ellipse(c(-0.05, 0.09), shape=A, radius=1.44, segments=1e3, col="red", lty=2,add=FALSE)
用它创建一个 data.table 并计算每个点到中心的距离 (point_x - center_x)² + (point_y - center_y )²:
dt <- data.table(tmp)
dt[,dist:={dx=x-center[1];dy=y-center[2];dx*dx+dy*dy}]
按距离对顶点排序：
setorder(dt,dist)
获取最小和最大点数：
> tail(dt,2)
x y dist
1: 4.990415 -6.138039 64.29517
2: -5.110415 6.318039 64.29517
> head(dt,2)
x y dist
1: 4.045722 3.41267 27.89709
2: -4.165722 -3.23267 27.89709
不要添加太多段，否则前两个值将是非常接近而不是相反的两个点。
从视觉效果来看，这听起来不太准确：
出于实用目的，@Tensibai 的回答可能已经足够好了。只需为 segments 参数使用足够大的值，以便这些点可以很好地近似于真实顶点。
如果你想要更严格的东西，你可以求解沿椭圆的位置，maximises/minimises 到中心的距离，由角度参数化。由于形状矩阵的存在，这比仅采用 angle={0, pi/2, pi, 3pi/2} 更复杂。不过也不算太难:
# location along the ellipse
# linear algebra lifted from the code for ellipse()
ellipse.loc <- function(theta, center, shape, radius)
{
vert <- cbind(cos(theta), sin(theta))
Q <- chol(shape, pivot=TRUE)
ord <- order(attr(Q, "pivot"))
t(center + radius*t(vert %*% Q[, ord]))
}
# distance from this location on the ellipse to the center
ellipse.rad <- function(theta, center, shape, radius)
{
loc <- ellipse.loc(theta, center, shape, radius)
(loc[,1] - center[1])^2 + (loc[,2] - center[2])^2
}
# ellipse parameters
center <- c(-0.05, 0.09)
A <- matrix(c(20.43, -8.59, -8.59, 24.03), nrow=2)
radius <- 1.44
# solve for the maximum distance in one hemisphere (hemi-ellipse?)
t1 <- optimize(ellipse.rad, c(0, pi - 1e-5), center=center, shape=A, radius=radius, maximum=TRUE)$m
l1 <- ellipse.loc(t1, center, A, radius)
# solve for the minimum distance
t2 <- optimize(ellipse.rad, c(0, pi - 1e-5), center=center, shape=A, radius=radius)$m
l2 <- ellipse.loc(t2, center, A, radius)
# other points obtained by symmetry
t3 <- pi + t1
l3 <- ellipse.loc(t3, center, A, radius)
t4 <- pi + t2
l4 <- ellipse.loc(t4, center, A, radius)
# plot everything
MASS::eqscplot(center[1], center[2], xlim=c(-7, 7), ylim=c(-7, 7), xlab="", ylab="")
ellipse(center, A, radius, col="red", lty=2)
points(rbind(l1, l2, l3, l4), cex=2, col="blue", lwd=2)
我知道这个问题已经解决了，但实际上有一个超级优雅的解决方案，只需几行如下。这种计算是精确的，没有任何形式的数值优化。
## target covariance matrix
A <- matrix(c(20.43, -8.59,-8.59, 24.03), nrow = 2)
E <- eigen(A, symmetric = TRUE) ## symmetric eigen decomposition
U <- E[[2]] ## eigen vectors, i.e., rotation matrix
D <- sqrt(E[[1]]) ## root eigen values, i.e., scaling factor
r <- 1.44 ## radius of original circle
Z <- rbind(c(r, 0), c(0, r), c(-r, 0), c(0, -r)) ## original vertices on major / minor axes
Z <- tcrossprod(Z * rep(D, each = 4), U) ## transformed vertices on major / minor axes
# [,1] [,2]
#[1,] -5.055136 6.224212
#[2,] -4.099908 -3.329834
#[3,] 5.055136 -6.224212
#[4,] 4.099908 3.329834
C0 <- c(-0.05, 0.09) ## new centre
Z <- Z + rep(C0, each = 4) ## shift to new centre
# [,1] [,2]
#[1,] -5.105136 6.314212
#[2,] -4.149908 -3.239834
#[3,] 5.005136 -6.134212
#[4,] 4.049908 3.419834
为了解释背后的数学原理，我将分3步：
- 这个椭圆是从哪里来的？
- Cholesky 分解方法及其缺点。
- 特征分解方法及其自然解释。
这个椭圆是哪里来的？
在实践中，这个椭圆可以通过对单位圆的一些线性变换得到x ^ 2 + y ^ 2 = 1。
Cholesky分解法及其缺点
## initial circle
r <- 1.44
theta <- seq(0, 2 * pi, by = 0.01 * pi)
X <- r * cbind(cos(theta), sin(theta))
## target covariance matrix
A <- matrix(c(20.43, -8.59,-8.59, 24.03), nrow = 2)
R <- chol(A) ## Cholesky decomposition
X1 <- X %*% R ## linear transformation
Z <- rbind(c(r, 0), c(0, r), c(-r, 0), c(0, -r)) ## original vertices on major / minor axes
Z1 <- Z %*% R ## transformed coordinates
## different colour per quadrant
g <- floor(4 * (1:nrow(X) - 1) / nrow(X)) + 1
## draw ellipse
plot(X1, asp = 1, col = g)
points(Z1, cex = 1.5, pch = 21, bg = 5)
## draw circle
points(X, col = g, cex = 0.25)
points(Z, cex = 1.5, pch = 21, bg = 5)
## draw axes
abline(h = 0, lty = 3, col = "gray", lwd = 1.5)
abline(v = 0, lty = 3, col = "gray", lwd = 1.5)
我们看到线性变换矩阵 R 似乎没有自然解释。圆的原始顶点不映射到椭圆的顶点。
特征分解方法及其自然解释
## initial circle
r <- 1.44
theta <- seq(0, 2 * pi, by = 0.01 * pi)
X <- r * cbind(cos(theta), sin(theta))
## target covariance matrix
A <- matrix(c(20.43, -8.59,-8.59, 24.03), nrow = 2)
E <- eigen(A, symmetric = TRUE) ## symmetric eigen decomposition
U <- E[[2]] ## eigen vectors, i.e., rotation matrix
D <- sqrt(E[[1]]) ## root eigen values, i.e., scaling factor
r <- 1.44 ## radius of original circle
Z <- rbind(c(r, 0), c(0, r), c(-r, 0), c(0, -r)) ## original vertices on major / minor axes
## step 1: re-scaling
X1 <- X * rep(D, each = nrow(X)) ## anisotropic expansion to get an axes-aligned ellipse
Z1 <- Z * rep(D, each = 4L) ## vertices on axes
## step 2: rotation
Z2 <- tcrossprod(Z1, U) ## rotated vertices on major / minor axes
X2 <- tcrossprod(X1, U) ## rotated ellipse
## different colour per quadrant
g <- floor(4 * (1:nrow(X) - 1) / nrow(X)) + 1
## draw rotated ellipse and vertices
plot(X2, asp = 1, col = g)
points(Z2, cex = 1.5, pch = 21, bg = 5)
## draw axes-aligned ellipse and vertices
points(X1, col = g)
points(Z1, cex = 1.5, pch = 21, bg = 5)
## draw original circle
points(X, col = g, cex = 0.25)
points(Z, cex = 1.5, pch = 21, bg = 5)
## draw axes
abline(h = 0, lty = 3, col = "gray", lwd = 1.5)
abline(v = 0, lty = 3, col = "gray", lwd = 1.5)
## draw major / minor axes
segments(Z2[1,1], Z2[1,2], Z2[3,1], Z2[3,2], lty = 2, col = "gray", lwd = 1.5)
segments(Z2[2,1], Z2[2,2], Z2[4,1], Z2[4,2], lty = 2, col = "gray", lwd = 1.5)
在这里我们看到，在变换的两个阶段，顶点仍然映射到顶点。正是基于这样 属性 我们在一开始就给出了巧妙的解决方案。