Добавление значений самозапуска к регрессии nls в R - PullRequest
Новая
       51

Добавление значений самозапуска к регрессии nls в R

1 голос
/ 16 января 2020

У меня есть существующий код для подбора сигмоидной кривой к данным в R. Как я могу использовать selfstart (или другой метод) для автоматического поиска начальных значений для регрессии? 

sigmoid = function(params, x) {
  params[1] / (1 + exp(-params[2] * (x - params[3])))
}

dataset = data.frame("x" = 1:53, "y" =c(0,0,0,0,0,0,0,0,0,0,0,0,0,0.1,0.18,0.18,0.18,0.33,0.33,0.33,0.33,0.41,0.41,0.41,0.41,0.41,0.41,0.5,0.5,0.5,0.5,0.68,0.58,0.58,0.68,0.83,0.83,0.83,0.74,0.74,0.74,0.83,0.83,0.9,0.9,0.9,1,1,1,1,1,1,1) )

x = dataset$x
y = dataset$y

# fitting code
fitmodel <- nls(y~a/(1 + exp(-b * (x-c))), start=list(a=1,b=.5,c=25))

# visualization code
# get the coefficients using the coef function
params=coef(fitmodel)

y2 <- sigmoid(params,x)
plot(y2,type="l")
points(y)

1 Ответ

3 голосов
/ 16 января 2020

Это распространенная (и интересная) проблема нелинейного подбора кривой.

Фон

Мы можем найти разумные начальные значения, если более внимательно посмотреть на функцию sigmoid

enter image description here

Сначала отметим, что

enter image description here

То же для больших значений x, функция приближается к a. Другими словами, в качестве начального значения для a мы можем выбрать значение y для наибольшего значения x. На языке R это означает y[which.max(x)].

Теперь, когда у нас есть начальное значение для a, нам нужно определиться с начальными значениями для b и c. Для этого мы можем использовать геометрию c серии

enter image description here

и расширить f(x) = y, сохранив только первые два члена

enter image description here

Теперь мы установим a = 1 (наше начальное значение для a), перестроим уравнение и возьмем логарифм с обеих сторон

enter image description here

Теперь мы можем подобрать линейную модель вида log(1 - y) ~ x для получения оценок наклона и смещения, которые, в свою очередь, дают начальные значения для b и c.

Реализация R

Давайте определим функцию, которая принимает в качестве аргумента значения x и y и возвращает list начальных значений параметров 

start_val_sigmoid <- function(x, y) {
    fit <- lm(log(y[which.max(x)] - y + 1e-6) ~ x)
    list(
        a = y[which.max(x)],
        b = unname(-coef(fit)[2]),
        c = unname(-coef(fit)[1] / coef(fit)[2]))
}

На основании данных для x и y, которые вы даете, мы получаем следующие начальные значения 

start_val_sigmoid(x, y)
#$a
#[1] 1
#
#$b
#[1] 0.2027444
#
#$c
#[1] 15.01613

Поскольку start_val_sigmoid возвращает list, мы можем использовать его выводить напрямую как аргумент start в nls

nls(y ~ a / ( 1 + exp(-b * (x - c))), start = start_val_sigmoid(x, y))
#Nonlinear regression model
#  model: y ~ a/(1 + exp(-b * (x - c)))
#   data: parent.frame()
#      a       b       c
# 1.0395  0.1254 29.1725
# residual sum-of-squares: 0.2119
#
#Number of iterations to convergence: 9
#Achieved convergence tolerance: 9.373e-06

Пример данных 

dataset = data.frame("x" = 1:53, "y" =c(0,0,0,0,0,0,0,0,0,0,0,0,0,0.1,0.18,0.18,0.18,0.33,0.33,0.33,0.33,0.41,0.41,0.41,0.41,0.41,0.41,0.5,0.5,0.5,0.5,0.68,0.58,0.58,0.68,0.83,0.83,0.83,0.74,0.74,0.74,0.83,0.83,0.9,0.9,0.9,1,1,1,1,1,1,1) )

x = dataset$x
y = dataset$y
...