2017-09-08 62 views
2

如何自動提取R^2對於整條曲線不理想的曲線的擬合線性部分?如何查找曲線的線性部分

例如 我有什麼:

data.lm

x y 
1 1 1 
2 2 8 
3 3 3 
4 4 4 
5 5 5 
6 6 6 
7 7 7 
8 8 5 
9 9 2 
10 10 7 

rg.lm < - 流明(Y〜X,data.lm) rg.lm

Coefficients: 
(Intercept)   x 
    3.7333  0.1939 

摘要(rg.lm)

Residuals: 
    Min  1Q Median  3Q  Max 
-3.4788 -1.1136 0.0061 1.2712 3.8788 

Coefficients: 
      Estimate Std. Error t value Pr(>|t|) 
(Intercept) 3.7333  1.6111 2.317 0.0491 * 
x    0.1939  0.2597 0.747 0.4765 

Residual standard error: 2.358 on 8 degrees of freedom 
Multiple R-squared: 0.06519, Adjusted R-squared: -0.05166 
F-statistic: 0.5579 on 1 and 8 DF, p-value: 0.4765 

我期待什麼:

data.lm.ex < - unknown.function(data.lm) data.lm.前

x y 
1 3 3 
2 4 4 
3 5 5 
4 6 6 
7 7 7 

Anoth呃例子來自真實數據:

data.lm

time OD 
1  0 2.175 
2 30 2.134 
3 60 2.189 
4 90 2.141 
5 120 2.854 
6 150 3.331 
7 180 3.642 
8 210 4.333 
9 240 4.987 
10 270 5.093 
11 300 4.943 
12 330 5.198 
13 360 4.804 

摘要(LM(data.lm))$ r.squared

[1] 0.8981063 

summary(lm(data.lm [4:9,]))$ r.squared

[1] 0.9886727 

正如上面所示,線4之間的間隔到9具有比整個曲線絕對更高R^2。你能否讓我知道用自動的方法來找到最高r^2的區間,並且至少有一定數量的點(由於2個點始終存在r^2 = 1.0)?

回答

2

這應該工作:

a <- cbind(1:10, c(1,8,3:7,5,2,7)) 
tmp <- rle(diff(a[,2])) 
ml <- max(tmp$lengths) 
i1 <- which(ml==tmp$lengths)[1] 

a[seq(i1,i1+ml),] 

更新

a <- data.frame(x=c(0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360), 
       y=c(2.175, 2.134, 2.189, 2.141, 2.854, 3.331, 3.642, 4.333, 4.987, 5.093, 4.943, 5.198, 4.804)) 

b <- diff(a[,2])/diff(a[,1]) 
b.k <- kmeans(b,3) 
b.max <- max(abs(b.k$centers)) 
b.v <- which(b.k$cluster == match(b.max, b.k$centers)) 

RES <- a[b.v,] 
plot(a) 
points(RES,pch=15) 
abline(coef(lm(y~x,RES)), col="red") 

enter image description here

改良版本:

library(zoo) 
a <- data.frame(x=c(0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360), 
       y=c(2.175, 2.134, 2.189, 2.141, 2.854, 3.331, 3.642, 4.333, 4.987, 5.093, 4.943, 5.198, 4.804)) 
f <- function (d) { 
    m <- lm(y~x, as.data.frame(d)) 
    return(coef(m)[2]) 
} 
co <- rollapply(a, 3, f, by.column=F) 
co.cl <- kmeans(co, 2) 
b.points <- which(co.cl$cluster == match(max(co.cl$centers), co.cl$centers))+1 
RES <- a[b.points,] 
plot(a) 
points(RES,pch=15,col="red") 
abline(lm(y~x,RES),col="blue") 

[an improved version]

+0

非常感謝。遺憾的是,這種解決方案僅適用於完全線性的配件。對於實際數據來說,R^2不可能是1.0。 –

+0

在編輯的問題中添加了一個真實的例子。在你方便的時候,你能幫我找出答案嗎? –

+0

當我看到您提供的這些新數據時,我會看到三個線性部分。你對哪一個感興趣? – akond