本文介绍了在R中使用nls()进行分段函数拟合的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我正在尝试将数据分为两部分.

以下是一些示例数据:

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) )

但是我得到了错误:

我从找到一条曲线来匹配数据

有什么想法我要去哪里吗?

解决方案

我没有一个很好的答案,但是我有一个一个答案.

(请参阅下面的编辑,以获取更高级的答案)

如果Cx足够小以致没有数据点适合ACy,或者Cx足够大以致没有数据点适合BCy到,QR分解矩阵将是奇异的,因为分别有CxACyCxBCy的许多不同值可以很好地拟合数据./p>

我通过防止安装Cx对此进行了测试.如果我在(说)Cx = mean(x)处修复了Cx,则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()解决了固定Cx时的问题,因此我现在可以找到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

对于f的最佳值,AB的值与yayb的值之间没有太大差异,但是存在一些差异.

(编辑-优雅答案)

将问题分为两个步骤,因此不再需要使用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
}

...给出:

I am trying to fit a two-part line to data.

Here's some sample data:

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)

I am hoping to be able to define the best fitting two part line (hand drawn example shown)

I then define a piecewise function that should find a two part linear function. The definition is based on the gradients of the two lines and their intercept with each other, which should completely define the lines.

# 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) )

However I get the error:

I got the basic method from Finding a curve to match data

Any ideas where I am going wrong?

解决方案

I don't have an elegant answer, but I do have an answer.

(SEE THE EDIT BELOW FOR A MORE ELEGANT ANSWER)

If Cx is small enough that there are no data points to fit A and Cy to, or if Cx is big enough that there are no data points to fit B and Cy to, the QR decomposition matrix will be singular because there will be many different values of Cx, A and Cy or Cx, B and Cy respectively that will fit the data equally well.

I tested this by preventing Cx from being fitted. If I fix Cx at (say) Cx = mean(x), nls() solves the problem without difficulty:

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))

... gives:

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

That led me to think that if I transformed Cx so that it could never go outside the range [min(x),max(x)], that might solve the problem. In fact, I'd want there to be at least three data points available to fit each of the "A" line and the "B" line, so Cx has to be between the third lowest and the third highest values of x. Using the atan() function with the appropriate arithmetic let me map a range [-inf,+inf] onto [0,1], so I got the code:

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))

Unfortunately, however, I still get the singular gradient matrix at initial parameters error from this code, so the problem is still over-parameterised. As @Henrik has suggested, the difference between the bilinear and single linear fit is not great for these data.

I can nevertheless get an answer for the bilinear fit, however. Since nls() solves the problem when Cx is fixed, I can now find the value of Cx that minimises the residual standard error by simply doing a one-dimensional minimisation using optimize(). Not a particularly elegant solution, but better than nothing:

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))

... gives output of:

[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

There isn't a huge difference between the values of A and B and ya and yb for the optimum value of f, but there is some difference.

(EDIT -- ELEGANT ANSWER)

Having separated the problem into two steps, it isn't necessary to use nls() any more. lm() works fine, as follows:

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
}

... which gives:

这篇关于在R中使用nls()进行分段函数拟合的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!

09-05 19:57
查看更多