Двойной линейный фитинг в R - PullRequest
Новая
       9

Двойной линейный фитинг в R

0 голосов
/ 24 мая 2019

У меня есть следующие данные: 

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)

Создано в 2019-05-24 пакетом Представить (v0.3.0)

Эти данные могут быть снабжены функцией двухлинейная , например: 

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.

Я знаю, что этого можно достичь с помощью *Пакет 1028 *, подобный следующему: 

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

Создан в 2019-05-24 с помощью пакета Представлять (v0.3.0)

Тем не менее, я хотел бы перевести эту функцию: 

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

Для достижения тех же результатов.Есть идеи?

1 Ответ

1 голос
/ 24 мая 2019

Я бы подошел, используя 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))

Создано в 2019-05-24 с помощью представительного пакета (v0.2.1)

PS Параметры, рассчитанные по подгонке, следующие: 

# 0.50036077 0.07611683 0.40440741 8.07065399

, которые находятся в тесном согласии с пакетом segmented: 

# 0.4970     0.0768     0.4032     8.07
...