7

我正在尝试将两部分的线拟合到数据中。

以下是一些示例数据：

x<-c(0.00101959664756622, 0.001929220749155, 0.00165657261751726, 
0.00182514724375389, 0.00161532360585458, 0.00126991061099209, 
0.00149545009309177, 0.000816386510029308, 0.00164402569283353, 
0.00128029006251656, 0.00206892841921455, 0.00132378793976235, 
0.000953143467154676, 0.00272964503695939, 0.00169743839571702, 
0.00286411493120396, 0.0016464862337286, 0.00155672067449593, 
0.000878271561566836, 0.00195872573138819, 0.00255412836538339, 
0.00126212428137799, 0.00106206607962734, 0.00169140916371657, 
0.000858015581562961, 0.00191955159274793, 0.00243104345247067, 
0.000871042201994687, 0.00229814264111745, 0.00226756341241083)

y<-c(1.31893118849162, 0.105150790530179, 0.412732029152914, 0.25589805483046, 
0.467147868109498, 0.983984462069833, 0.640007862668818, 1.51429617241365, 
0.439777145282391, 0.925550163462951, -0.0555942758921906, 0.870117027565708, 
1.38032147826294, -0.96757052387814, 0.346370836378525, -1.08032147826294, 
0.426215616848312, 0.55151485221263, 1.41306889485598, 0.0803478641720901, 
-0.86654892295057, 1.00422341998656, 1.26214517662281, 0.359512373951839, 
1.4835398594013, 0.154967053938309, -0.680501679226447, 1.44740598234453, 
-0.512732029152914, -0.359512373951839)

我希望能够定义最合适的两部分线（显示手绘示例）

阴谋

然后我定义了一个分段函数，它应该找到一个两部分的线性函数。该定义基于两条线的梯度及其相互截距，这应该完全定义线。

# A=gradient of first line segment
# B=gradient of second line segment
# Cx=inflection point x coord
# Cy=inflexion point y coord 

out_model <- nls(y ~ I(x <= Cx)*Cy-A*(Cx-x)+I(x > Cx)*Cy+B*(x), 
                  data = data.frame(x,y), 
                  start = c(A=-500,B=-500,Cx=0.0001,Cy=-1.5) )

但是我得到了错误：

nls(y ~ I(x <= Cx) * Cy - A * (Cx - x) + I(x > Cx) * Cy + B * 中的错误：奇异梯度

我从Finding a curve to match data中得到了基本方法

有什么想法我哪里出错了吗？

4

4 回答 4

4

我没有优雅的答案，但我有一个答案。

（请参阅下面的编辑以获得更优雅的答案）

如果Cx足够小以至于没有数据点可以拟合A和Cy到，或者如果Cx足够大以至于没有数据点可以拟合B和到，则 QR 分解矩阵将是奇异的，因为,和或Cy会有许多不同的值，并且分别可以很好地拟合数据。CxACyCxBCy

我通过防止Cx安装来测试这一点。如果我修复Cx(say) Cx = mean(x)，nls()可以毫无困难地解决问题：

nls(y ~ ifelse(x < mean(x),ya+A*x,yb+B*x), 
               data = data.frame(x,y), 
               start = c(A=-1000,B=-1000,ya=3,yb=0))

...给出：

Nonlinear regression model
  model:  y ~ ifelse(x < mean(x), ya + A * x, yb + B * x) 
   data:  data.frame(x, y) 
        A         B        ya        yb 
-1325.537 -1335.918     2.628     2.652 
 residual sum-of-squares: 0.06614

Number of iterations to convergence: 1 
Achieved convergence tolerance: 2.294e-08

这使我认为，如果我进行改造Cx以使其永远不会超出范围[min(x),max(x)]，那可能会解决问题。事实上，我希望至少有三个数据点可用于拟合“A”线和“B”线中的每一个，因此 Cx 必须介于 的第三低和第三高值之间x。使用atan()具有适当算术的函数让我将范围映射[-inf,+inf]到 上[0,1]，所以我得到了代码：

trans <- function(x) 0.5+atan(x)/pi
xs <- sort(x)
xlo <- xs[3]
xhi <- xs[length(xs)-2]
nls(y ~ ifelse(x < xlo+(xhi-xlo)*trans(f),ya+A*x,yb+B*x), 
               data = data.frame(x,y), 
               start = c(A=-1000,B=-1000,ya=3,yb=0,f=0))

然而不幸的是，我仍然singular gradient matrix at initial parameters从这段代码中得到错误，所以问题仍然是过度参数化的。正如@Henrik 所建议的那样，双线性和单线性拟合之间的差异对于这些数据来说并不是很大。

然而，我仍然可以得到双线性拟合的答案。由于nls()解决了问题 whenCx是固定的，我现在可以Cx通过简单地使用 进行一维最小化来找到最小化残差标准误差的值optimize()。不是一个特别优雅的解决方案，但总比没有好：

xs <- sort(x)
xlo <- xs[3]
xhi <- xs[length(xs)-2]
nn <- function(f) nls(y ~ ifelse(x < xlo+(xhi-xlo)*f,ya+A*x,yb+B*x), 
               data = data.frame(x,y), 
               start = c(A=-1000,B=-1000,ya=3,yb=0))
ssr <- function(f) sum(residuals(nn(f))^2)
f = optimize(ssr,interval=c(0,1))
print (f$minimum)
print (nn(f$minimum))
summary(nn(f$minimum))

...给出以下输出：

[1] 0.8541683
Nonlinear regression model
  model:  y ~ ifelse(x < xlo + (xhi - xlo) * f, ya + A * x, yb + B * x) 
   data:  data.frame(x, y) 
        A         B        ya        yb 
-1317.215  -872.002     2.620     1.407 
 residual sum-of-squares: 0.0414

Number of iterations to convergence: 1 
Achieved convergence tolerance: 2.913e-08 

Formula: y ~ ifelse(x < xlo + (xhi - xlo) * f, ya + A * x, yb + B * x)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
A  -1.317e+03  1.792e+01 -73.493  < 2e-16 ***
B  -8.720e+02  1.207e+02  -7.222 1.14e-07 ***
ya  2.620e+00  2.791e-02  93.854  < 2e-16 ***
yb  1.407e+00  3.200e-01   4.399 0.000164 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.0399 on 26 degrees of freedom

Number of iterations to convergence: 1

的和的值A和的最佳值之间没有很大的差异，但是有一些差异。Byaybf

（编辑——优雅的答案）

把问题分成两步，就没有必要再用nls()了。lm()工作正常，如下：

function (x,y) 
{
    f <- function (Cx) 
        {
        lhs <- function(x) ifelse(x < Cx,Cx-x,0)
        rhs <- function(x) ifelse(x < Cx,0,x-Cx)
        fit <- lm(y ~ lhs(x) + rhs(x))
        c(summary(fit)$r.squared, 
            summary(fit)$coef[1], summary(fit)$coef[2],
            summary(fit)$coef[3])
        }

    r2 <- function(x) -(f(x)[1])

    res <- optimize(r2,interval=c(min(x),max(x)))
    res <- c(res$minimum,f(res$minimum))

    best_Cx <- res[1]
    coef1 <- res[3]
    coef2 <- res[4]
    coef3 <- res[5]
    plot(x,y)
    abline(coef1+best_Cx*coef2,-coef2) #lhs  
    abline(coef1-best_Cx*coef3,coef3)  #rs
}

... 这使：

在此处输入图像描述

于 2013-04-08T11:12:08.117 回答
3

如果断点已知，则可以使用线性回归

来自“使用 R 的实际回归和方差分析”的断棒回归

朱利安·J·法拉维

2000 年 12 月

k <- 0.0025

lhs <- function(x) ifelse(x < k,k-x,0)
rhs <- function(x) ifelse(x < k,0,x-k)
fit <- lm(y ~ lhs(x) + rhs(x))
于 2013-04-08T11:32:39.093 回答
3

该软件包segmented专为此类问题而设计。

首先，创建一个常规线性回归lm：

linmod <- lm(y ~ x)
summary(linmod)

这给了我们：

Call:
lm(formula = y ~ x)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.108783 -0.025432 -0.006484  0.040092  0.088638 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.630e+00  2.732e-02   96.28   <2e-16 ***
x           -1.326e+03  1.567e+01  -84.63   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04869 on 28 degrees of freedom
Multiple R-squared:  0.9961,    Adjusted R-squared:  0.996 
F-statistic:  7163 on 1 and 28 DF,  p-value: < 2.2e-16

接下来，我们使用线性模型来生成具有 1 个断点的分段模型：

segmod <- segmented(linmod, seg.Z = ~x)
summary(segmod)

分段模型提供了更好的 r 平方：

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = linmod, seg.Z = ~x)

Estimated Break-Point(s):
   Est. St.Err 
 0.003  0.000 

Meaningful coefficients of the linear terms:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.659e+00  2.882e-02  92.239   <2e-16 ***
x           -1.347e+03  1.756e+01 -76.742   <2e-16 ***
U1.x         5.167e+02  4.822e+02   1.072       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04582 on 26 degrees of freedom
Multiple R-Squared: 0.9968,  Adjusted R-squared: 0.9964 

Convergence attained in 3 iterations with relative change 0

您可以检查绘图、截距和斜率：

plot(segmod)
intercept(segmod)
slope(segmod)
于 2018-05-09T12:35:51.610 回答
1

感谢 Henrik 让我走上了正确的道路！这是一个更完整且相对优雅的解决方案，带有一个简单的情节：

range_x<-max(x)-min(x)
intervals=1000
coef1=c()
coef2=c()
coef3=c()
r2=c()

for (i in 1:intervals)  
{
Cx<-min(x)+(i-1)*(range_x/intervals)
lhs <- function(x) ifelse(x < Cx,Cx-x,0)
rhs <- function(x) ifelse(x < Cx,0,x-Cx)
fit <- lm(y ~ lhs(x) + rhs(x))
coef1[i]<-summary(fit)$coef[1]
coef2[i]<-summary(fit)$coef[2]
coef3[i]<-summary(fit)$coef[3]
r2[i]<-summary(fit)$r.squared
}
best_r2<-max(r2)                             # get best r squared
pos<-which.max(r2)                                          
best_Cx<-min(x)+(pos-1)*(range_x/intervals)  # get Cx for best r2

plot(x,y)
abline(coef1[pos]+best_Cx*coef2[pos],-coef2[pos]) #lhs  
abline(coef1[pos]-best_Cx*coef3[pos],coef3[pos])  #rs

在此处输入图像描述

于 2013-04-09T02:48:32.623 回答