r - 使用 nls() 进行非线性拟合在初始参数估计时给了我奇异的梯度矩阵。为什么？

Question

这是我第一次尝试在 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

任何帮助是极大的赞赏。

Answer 1

我假设你得到了一个长度为 len 的观察向量。如您的示例中绘制的那些，并且您希望识别 k跳跃和k跳跃大小。 (或者也许我误解了你；但你并没有真正说出你想要达到的目标。)
下面我将使用本地搜索绘制一个解决方案。我从您的示例数据开始:

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