Нелинейная аппроксимация с помощью nls () дает мне особую градиентную матрицу при начальных оценках параметров. Зачем? - PullRequest
3 голосов
/ 16 мая 2019

Это моя первая попытка подобрать нелинейную модель в R, поэтому, пожалуйста, потерпите меня.

Проблема

Я пытаюсь понять, почему nls() дает мне этоошибка:

Error in nlsModel(formula, mf, start, wts): singular gradient matrix at initial parameter estimates

Гипотезы

Из того, что я прочитал из других вопросов здесь, в SO, это может быть либо потому, что:

  • моя модель прерывистая,или
  • моя модель переопределена, или
  • неправильный выбор начальных значений параметров

Поэтому я обращаюсь за помощью о том, как преодолеть эту ошибку.Могу ли я изменить модель и все еще использовать nls(), или мне нужно использовать nls.lm из пакета minpack.lm, как я читал в другом месте?

Мой подход

Вотнекоторые подробности о модели:

  • модель представляет собой прерывистую функцию, тип функции лестница (см. график ниже)
  • в целом,количество шагов в модели может быть переменным, но они фиксированы для конкретного события подгонки

MWE, которое показывает проблему

Краткое описание кода MWE

  • step_fn(x, min = 0, max = 1): функция, которая возвращает 1 в пределах интервала (min, max] и 0 в противном случае; извините за имя, теперь я понимаю, что это не совсемШаговая функция ... interval_fn() было бы более подходящим, я думаю.
  • staircase(x, dx, dy): сумма step_fn() функций. dx - вектор ширины для шагов , т.е. max - min, а dy - это приращение y для каждого шага .
  • staircase_formula(n = 1L): генерирует formulaобъект, представляющий модель, моделируемую функцией staircase() (для использования с функцией nls()).
  • обратите внимание, что я использую пакеты purrr и glue в приведенном ниже примере.

Код

step_fn <- function(x, min = 0, max = 1) {

  y <- x
  y[x > min & x <= max] <- 1
  y[x <= min] <- 0
  y[x > max] <- 0

  return(y)
}

staircase <- function(x, dx, dy) {

  max <- cumsum(dx)
  min <- c(0, max[1:(length(dx)-1)])
  step <- cumsum(dy)

  purrr::reduce(purrr::pmap(list(min, max, step), ~ ..3 * step_fn(x, min = ..1, max = ..2)), `+`)
}


staircase_formula <- function(n = 1L) {

  i <- seq_len(n)
  dx <- sprintf("dx%d", i)

  min <-
    c('0', purrr::accumulate(dx[-n], .f = ~ paste(.x, .y, sep = " + ")))
  max <- purrr::accumulate(dx, .f = ~ paste(.x, .y, sep = " + "))

  lhs <- "y"
  rhs <-
    paste(glue::glue('dy{i} * step_fn(x, min = {min}, max = {max})'),
          collapse  = " + ")

  sc_form <- as.formula(glue::glue("{lhs} ~ {rhs}")) 

  return(sc_form)
}


x <- seq(0, 10, by = 0.01)
y <- staircase(x, c(1,2,2,5), c(2,5,2,1)) + rnorm(length(x), mean = 0, sd = 0.2)

plot(x = x, y = y)
lines(x = x, y = staircase(x, dx = c(1,2,2,5), dy = c(2,5,2,1)), col="red")


my_data <- data.frame(x = x, y = y)
my_model <- staircase_formula(4)
params <- list(dx1 = 1, dx2 = 2, dx3 = 2, dx4 = 5,
               dy1 = 2, dy2 = 5, dy3 = 2, dy4 = 1)

m <- nls(formula = my_model, start = params, data = my_data)
#> Error in nlsModel(formula, mf, start, wts): singular gradient matrix at initial parameter estimates

Любая помощь очень ценится.

Ответы [ 2 ]

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

Я предполагаю, что вам дан вектор наблюдений длины len, как те, что изображены в вашем примере, и вы хотите определить k скачков и k размеров скачков.(Или, может быть, я вас неправильно понял, но вы на самом деле не сказали, чего хотите достичь.) Ниже я нарисую решение, используя Local Search.Я начну с данных вашего примера:

x <- seq(0, 10, by = 0.01)
y <- staircase(x,
               c(1,2,2,5),
               c(2,5,2,1)) + rnorm(length(x), mean = 0, sd = 0.2)

Решением является список позиций и размеров прыжков.Обратите внимание, что я использую векторы для хранения этих данных, так как определение переменных становится громоздким, скажем, при 20 скачках.

Пример (случайное) решение:

k <- 5   ## number of jumps
len <- length(x)

sol <- list(position = sample(len, size = k),
            size = runif(k))

## $position
## [1]  89 236 859 885 730
## 
## $size
## [1] 0.2377453 0.2108495 0.3404345 0.4626004 0.6944078

Нам нужноцелевая функция для расчета качества решения.Я также определяю простую вспомогательную функцию stairs, которая используется целевой функцией.Целевая функция abs_diff вычисляет среднюю абсолютную разницу между подобранными рядами (как определено решением) и y.

stairs <- function(len, position, size) {
    ans <- numeric(len)
    ans[position] <- size
    cumsum(ans)
}

abs_diff <- function(sol, y, stairs, ...) {
    yy <- stairs(length(y), sol$position, sol$size)
    sum(abs(y - yy))/length(y)
}

Теперь идет ключевой компонент для локального поиска: функция соседства, котораяиспользуется для развития решения.Функция соседства принимает решение и слегка его меняет.Здесь он либо выберет позицию или размер и немного его изменит.

neighbour <- function(sol, len, ...) {
    p <- sol$position
    s <- sol$size

    if (runif(1) > 0.5) {
        ## either move one of the positions ...
        i <- sample.int(length(p),  size = 1)
        p[i] <- p[i] + sample(-25:25, size = 1)
        p[i] <- min(max(1, p[i]), len)        
    } else {
        ## ... or change a jump size
        i <- sample.int(length(s), size = 1)
        s[i] <- s[i] + runif(1, min = -s[i], max = 1)
    }

    list(position = p, size = s)
}

Пример вызова: здесь новое решение имеет свой первый размер прыжкаизменено.

## > sol
## $position
## [1]  89 236 859 885 730
## 
## $size
## [1] 0.2377453 0.2108495 0.3404345 0.4626004 0.6944078
## 
## > neighbour(sol, len)
## $position
## [1]  89 236 859 885 730
## 
## $size
## [1] 0.2127044 0.2108495 0.3404345 0.4626004 0.6944078

Мне осталось запустить локальный поиск.

library("NMOF")
sol.ls <- LSopt(abs_diff,
                list(x0 = sol, nI = 50000, neighbour = neighbour),
                stairs = stairs,
                len = len,
                y = y)

Мы можем построить решение: выделенная линия показана синим.

plot(x, y)
lines(x, stairs(len, sol.ls$xbest$position, sol.ls$xbest$size),
      col = "blue", type = "S")

Data and fitted line (in blue)

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

Попробуйте вместо DE:

library(NMOF)
 yf= function(params,x){
   dx1 = params[1]; dx2 = params[2]; dx3 = params[3]; dx4 = params[4];
   dy1 = params[5]; dy2 = params[6]; dy3 = params[7]; dy4 = params[8]
   dy1 * step_fn(x, min = 0, max = dx1) + dy2 * step_fn(x, min = dx1, 
               max = dx1 + dx2) + dy3 * step_fn(x, min = dx1 + dx2, max = dx1 + 
               dx2 + dx3) + dy4 * step_fn(x, min = dx1 + dx2 + dx3, max = dx1 + 
               dx2 + dx3 + dx4)
 }

 algo1 <- list(printBar = FALSE,
               nP  = 200L,
               nG  = 1000L,
               F   = 0.50,
               CR  = 0.99,
               min = c(0,1,1,4,1,4,1,0),
               max = c(2,3,3,6,3,6,3,2))

 OF2 <- function(Param, data) { #Param=paramsj data=data2
   x <- data$x
   y <- data$y
   ye <- data$model(Param,x)
   aux <- y - ye; aux <- sum(aux^2)
   if (is.na(aux)) aux <- 1e10
   aux
 }

 data5 <- list(x = x, y = y,  model = yf, ww = 1)
 system.time(sol5 <- DEopt(OF = OF2, algo = algo1, data = data5))
 sol5$xbest
 OF2(sol5$xbest,data5)

 plot(x,y)
 lines(data5$x,data5$model(sol5$xbest, data5$x),col=7,lwd=2)

#>  sol5$xbest
#[1]   1.106396  12.719182  -9.574088  18.017527   3.366852   8.721374 -19.879474   1.090023
#>  OF2(sol5$xbest,data5)
#[1] 1000.424

enter image description here

...