在R中对3D点云拟合一条直线

Question

在R中对3D点云拟合一条直线

3

我有一些点云数据（xyz坐标），需要拟合成线性模型。我打算使用lm()函数进行拟合。

这是我的尝试过程：

library(scatterplot3d)

# example points
x <- c(1,4,3,6,2,5)
y <- c(2,2,4,3,5,9)
z <- c(1,3,5,9,2,2)

# plot    
s <- scatterplot3d(x,y,z, type="b")

# fit the model
ff = lm(z ~ x + y) ## in ff$coefficients are the line paramters z, mx, ny

# create coordinates for a short line (the fit) to plot
llx = c(min(x), max(x))
lly = c(min(y), max(y))
llz = c(
  ff$coefficients[[1]] + llx[1] * ff$coefficients[[2]] + lly[1] * ff$coefficients[[3]],
  ff$coefficients[[1]] + llx[2] * ff$coefficients[[2]] + lly[2] * ff$coefficients[[3]]
)

## create 2d coordinates to place in scatterplot
p0 <- s$xyz.convert(llx[1],lly[1],llz[1])
p1 <- s$xyz.convert(llx[2],lly[2],llz[2])

# draw line
segments(p0$x,p0$y,p1$x,p1$y,lwd=2,col=2)

尽管红线看起来像是很合适的拟合，但我并不确定它是否真的合适。如果你旋转图表，它看起来并不是很好。

for(i in seq(from=30, to=60, by=1)){
  s <- scatterplot3d(x,y,z, type="b", angle=i)
  segments(p0$x,p0$y,p1$x,p1$y,lwd=2,col=2)  
  Sys.sleep(0.1)
}

这是否只是由于线的二维投影？您能否以某种方式更新坐标？我尝试给 $xyz.convert() 函数添加 "angle" 属性，但没有成功。

另外，当我只使用两个示例点时，拟合失败。

x <- c(1,4)
y <- c(2,5)
z <- c(1,3)

我希望确认我是否正确使用lm()函数。谢谢！

[编辑]

我了解到，lm()函数基于我提供的模型（z~x+y）将一个平面拟合到数据上。这不是我想要的。实际上，我完全误解了lm()函数。即使对于二维数据也是如此。例如，lm(y~x)试图最小化拟合与数据之间的垂直距离。但是，我希望将数据视为完全独立的（空间数据），并最小化拟合与数据之间的垂线距离（参见这里的第一段：http://mathpages.com/home/kmath110.htm）。

标记为正确答案的回答正是这样做的。这个原则被称为“主成分分析”。

- agoldev

回归使用是正确的。问题在于绘制回归线。llx和lly的有序对可能并不总是对应于输入。就像上面的有序对llx和lly是(1,2)和(6,9)。虽然你通过运气得到了(1,2)对，但是(6,9)对不是你的输入数据集的一部分。因此，你得到了不正确的回归线。解决方案是使用输入数据集中的点的有序对来绘制回归线。 - 9Heads

llx和lly只是示例点的最大范围。我认为我也可以使用llx = c（0,1000000000）; lly = llx。但我使用了一条较短的线。我在这里得到llx = c（1,6）和lly = c（2,9）。所以，我认为你的答案没有解决问题，对吗？ - agoldev

1

想象一条与输入模型对应的直线（而不是线段）。尽管范围可能从“-Inf”到“Inf”。该直线对应于特定的三元组（x，y，z）（在您的情况下，是模型的输入和其他类似的输入和输出）。您不能期望该直线沿着您可以使用x和y独立地进行并相应的z所可能产生的所有三元组。然后它将覆盖整个“2D平面”。我所说的是，“x”和“y”是相关的（对应于相同的观察）。 - 9Heads

说实话，我不知道模型（z〜x + y）代表什么。我以为它是类似于2D线条的y = mx + n的三维表达式。这就是所有麻烦的开始之处;)对于点云的线性拟合，什么将是模型？ - agoldev

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- cuttlefish44 · Accepted Answer

lm(z ~ x + y)的拟合点并不是一条直线，而是一个平面。你的线段确实属于这个平面。

s <- scatterplot3d(x,y,z, type="b")
s$plane3d(ff)
segments(p0$x,p0$y,p1$x,p1$y,lwd=2,col=2) 

# rgl
library(rgl)
plot3d(x, y, z, type="s", rad=0.1)
planes3d(ff$coef[2], ff$coef[3], -1, ff$coef[1], col = 4, alpha = 0.3)
segments3d(llx, lly, llz, lwd=2, col=2)

[编辑]

您想要将三维数据拟合成一条直线，换句话说，将三维数据概括为一维数据。我认为这条直线由主成分分析的第一组成部分组成（即平均值 + t * PC1，该直线最小化总二乘误差）。我参考了 "R邮件帮助：将三维数据点拟合成一条直线" 和 "MathWorks：使用主成分分析进行正交回归拟合"。

x <- c(1,4,3,6,2,5)
y <- c(2,2,4,3,5,9)
z <- c(1,3,5,9,2,2)

xyz <- data.frame(x = x, y = y, z = z)
N <- nrow(xyz) 

mean_xyz <- apply(xyz, 2, mean)
xyz_pca   <- princomp(xyz) 
dirVector <- xyz_pca$loadings[, 1]   # PC1

xyz_fit <- matrix(rep(mean_xyz, each = N), ncol=3) + xyz_pca$score[, 1] %*% t(dirVector) 

t_ends <- c(min(xyz_pca$score[,1]) - 0.2, max(xyz_pca$score[,1]) + 0.2)  # for both ends of line
endpts <- rbind(mean_xyz + t_ends[1]*dirVector, mean_xyz + t_ends[2]*dirVector)

library(scatterplot3d) 
s3d <- scatterplot3d(xyz, type="b")
s3d$points3d(endpts, type="l", col="blue", lwd=2)
for(i in 1:N) s3d$points3d(rbind(xyz[i,], xyz_fit[i,]), type="l", col="green3", lty=2)

library(rgl)
plot3d(xyz, type="s", rad=0.1)
abclines3d(mean_xyz, a = dirVector, col="blue", lwd=2)     # mean + t * direction_vector
for(i in 1:N) segments3d(rbind(xyz[i,], xyz_fit[i,]), col="green3")