【发布时间】:2011-10-30 04:36:33
【问题描述】:
我从 R 脚本中的 lm 调用中获得了两个线性拟合。比如……
fit1 <- lm(y1 ~ x1)
fit2 <- lm(y2 ~ x2)
我想找到这两条线(fit1 和 fit2)相交的 (x,y) 点(如果它们完全相交的话)。
【问题讨论】:
标签: r intersection lm
我从 R 脚本中的 lm 调用中获得了两个线性拟合。比如……
fit1 <- lm(y1 ~ x1)
fit2 <- lm(y2 ~ x2)
我想找到这两条线(fit1 和 fit2)相交的 (x,y) 点(如果它们完全相交的话)。
【问题讨论】:
标签: r intersection lm
我有点惊讶没有内置函数。
这是一个基本函数(用于 lm 结果),使用与上面的 Tommy 相同的通用方法。这使用“y=mx+b”形式的两条线的简单替换方法来找到 y 处的公共交点 (y1=y2 ; m1*x + b1 = m2*x + b2) 并求解 x:
函数定义
# Linear model Intercept function
lmIntx <- function(fit1, fit2, rnd=2) {
b1<- fit1$coefficient[1] #y-int for fit1
m1<- fit1$coefficient[2] #slope for fit1
b2<- fit2$coefficient[1] #y-int for fit2
m2<- fit2$coefficient[2] #slope for fit2
if(m1==m2 & b1==b2) {print("Lines are identical")
} else if(m1==m2 & b1 != b2) {print("Lines are parallel")
} else {
x <- (b2-b1)/(m1-m2) #solved general equation for x
y <- m1*x + b1 #plug in the result
data.frame(x=round(x, rnd), y=round(y, rnd))
}
}
测试:
line1 <- data.frame(x=c(0,1), y=c(0,2))
line2 <- data.frame(x=c(0,1), y=c(1,3))
line3 <- data.frame(x=c(0,1), y=c(1,5))
lmIntx(lm(line1$y~line1$x), lm(line2$y~line2$x))
[1] "Lines are parallel"
lmIntx(lm(line1$y~line1$x), lm(line1$y~line1$x))
[1] "Lines are identical"
lmIntx(lm(line1$y~line1$x), lm(line3$y~line3$x))
x y
(Intercept) -0.5 -1
【讨论】:
避免几何的一种方法是将方程重新参数化为:
y1 = m1 * (x1 - x0) + y0
y2 = m2 * (x2 - x0) + y0
根据它们的交点(x0, y0),然后使用nls 一次执行两者的拟合,以便x0 和y0 的返回值给出结果:
# test data
set.seed(123)
x1 <- 1:10
y1 <- -5 + x1 + rnorm(10)
x2 <- 1:10
y2 <- 5 - x1 + rnorm(10)
g <- rep(1:2, each = 10) # first 10 are from x1,y1 and second 10 are from x2,y2
xx <- c(x1, x2)
yy <- c(y1, y2)
nls(yy ~ ifelse(g == 1, m1 * (xx - x0) + y0, m2 * (xx - x0) + y0),
start = c(m1 = -1, m2 = 1, y0 = 0, x0 = 0))
编辑:请注意,xx<-... 和 yy<-... 行是新的,nls 行已根据这些行指定并更正。
【讨论】:
algorithm = "port", lower = ...whatever..., upper = ...whatever... 按照?nls。
x0 介于2 和4 之间:nls(c(y1, y2) ~ ifelse(g == 1, b1 * (x1 - x0) + y0, b2 * (x2 - x0) + y0), start = c(b1 = -1, b2 = 1, y0 = 0, x0 = 3), algorithm = "port", lower = c(b1 = -Inf, b2 = -Inf, y0 = -Inf, x0 = 2), upper = c(b1 = Inf, b2 = Inf, y0 = Inf, x0 = 4))
这里有一些高中几何 ;-)
# First two models
df1 <- data.frame(x=1:50, y=1:50/2+rnorm(50)+10)
m1 <- lm(y~x, df1)
df2 <- data.frame(x=1:25, y=25:1*2+rnorm(25)-10)
m2 <- lm(y~x, df2)
# Plot them to show the intersection visually
plot(df1)
points(df2)
# Now calculate it!
a <- coef(m1)-coef(m2)
c(x=-a[[1]]/a[[2]], y=coef(m1)[[2]]*x + coef(m1)[[1]])
或者,通过@Dwin 简化基于solve 的解决方案:
cm <- rbind(coef(m1),coef(m2)) # Coefficient matrix
c(-solve(cbind(cm[,2],-1)) %*% cm[,1])
# [1] 12.68034 16.57181
【讨论】:
solve。虽然也许这更像是一个大学解决方案? ;)
solve 解决方案的更紧凑的变体。