R中的双线性拟合答案

【问题标题】：Double linear fitting in RR中的双线性拟合
【发布时间】：2019-05-24 10:10:47
【问题描述】：

我有以下数据：

df <- structure(list(x = c(0, 2.5, 5, 7.5, 10, 12.5, 15), 
                     y = c(0.51,0.71, 0.8, 1.12, 2.05, 3.23, 4.45)), 
                class = c("tbl_df", "tbl","data.frame"), row.names = c(NA, -7L))

df
#>      x    y
#> 1  0.0 0.51
#> 2  2.5 0.71
#> 3  5.0 0.80
#> 4  7.5 1.12
#> 5 10.0 2.05
#> 6 12.5 3.23
#> 7 15.0 4.45

plot(df)

^{由reprex package (v0.3.0) 于 2019 年 5 月 24 日创建}

这些数据可以用 双线性 函数拟合，如下所示：

if(x < bkp) {
  y <- i1 + s1 * x
} else {
  y <- (i1 + s1 * bkp) + s2 * (x - bkp)
}

其中bkp 是断点（介于7.5 和10 之间），i1 是y-intercept，s1 和s2 是slopes。

我知道这可以使用segmented 包来实现，如下所示：

library(segmented)

df <- structure(list(x = c(0, 2.5, 5, 7.5, 10, 12.5, 15), 
                     y = c(0.51,0.71, 0.8, 1.12, 2.05, 3.23, 4.45)), 
                class = c("tbl_df", "tbl","data.frame"), row.names = c(NA, -7L))

fit_df <- lm(y ~ x, data = df)

segmented(fit_df)
#> Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts
#> argument ignored
#> Call: segmented.lm(obj = fit_df)
#> 
#> Meaningful coefficients of the linear terms:
#> (Intercept)            x         U1.x  
#>      0.4970       0.0768       0.4032  
#> 
#> Estimated Break-Point(s):
#> psi1.x  
#>   8.07

^{由reprex package (v0.3.0) 于 2019 年 5 月 24 日创建}

但是，我想翻译这个函数：

if(x < bkp) {
  y <- i1 + s1 * x
} else {
  y <- (i1 + s1 * bkp) + s2 * (x - bkp)
}

为了达到同样的效果。有什么想法吗？

【问题讨论】：

如果它可能有一些用处，我可以很好地拟合一个简单的指数方程“y = a * exp(b * x)”，拟合参数 a = 4.0275246876414617E-01 和 b = 1.6137389391657100E-01

标签： r linear-regression

【解决方案1】：

我会使用 optim 和成本函数来拟合它。首先，我创建数据框。

# Data frame
df <- structure(list(x = c(0, 2.5, 5, 7.5, 10, 12.5, 15), 
                     y = c(0.51,0.71, 0.8, 1.12, 2.05, 3.23, 4.45)), 
                class = c("tbl_df", "tbl","data.frame"), row.names = c(NA, -7L))

接下来，我定义模型函数。请注意，我使用ifelse 简洁地将函数部分切换到断点右侧。

# Linear model with break point
model <- function(x, par){
  par[1] + par[2] * x + ifelse(x > par[4], par[3] * (x - par[4]), 0)
}

然后，我定义成本函数。这将计算平方残差的总和，并将其最小化以拟合模型。

# Cost function
cost <- function(par, df_data = df){
  sum((df_data$y - model(df_data$x, par))^2)
}

我致电optim 以最小化成本函数并绘制结果。

# Minimise cost function
fit <- optim(c(0, 0.1, 2, 7), cost)

# Plot results
plot(df)
lines(0:15, model(0:15, fit$par))

^{由reprex package (v0.2.1) 于 2019 年 5 月 24 日创建}

PS拟合估计的参数如下：

# 0.50036077 0.07611683 0.40440741 8.07065399

与segmented 包非常一致：

# 0.4970     0.0768     0.4032     8.07

【讨论】：