L'ajustement non linéaire avec nls() me donne une matrice de gradient singulière aux estimations initiales des paramètres. Pourquoi ?

Question

Il s'agit de ma première tentative d'ajustement d'un modèle non linéaire dans R, alors soyez indulgent avec moi.

Problème

J'essaie de comprendre pourquoi nls() me donne cette erreur :

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

Hypothèses

D'après ce que j'ai lu dans d'autres questions ici à SO, cela pourrait être soit parce que :

• mon modèle est discontinu, ou
• mon modèle est surdéterminé, ou
• mauvais choix des valeurs des paramètres de départ

J'appelle donc à l'aide pour savoir comment surmonter cette erreur. Puis-je changer le modèle et continuer à utiliser nls() ou dois-je utiliser nls.lm de la minpack.lm paquet, comme je l'ai lu ailleurs ?

Mon approche

Voici quelques détails sur le modèle :

• le modèle est une fonction discontinue, une sorte de escalier type de fonction (voir le graphique ci-dessous)
• en général, le nombre de étapes dans le modèle peuvent être variables, mais elles sont fixes pour un événement d'ajustement spécifique.

MWE qui montre le problème

Brève explication du code MWE

• step_fn(x, min = 0, max = 1) : fonction qui renvoie 1 dans l'intervalle ( min , max ] et 0 autrement ; désolé pour le nom, je réalise maintenant que ce n'est pas vraiment une fonction d'étape... interval_fn() serait plus approprié, je suppose.
• staircase(x, dx, dy) : une somme de step_fn() fonctions. dx est un vecteur de largeurs pour les étapes c'est-à-dire max - min y dy est l'augmentation de y pour chaque étape .
• staircase_formula(n = 1L) : génère un formula qui représente le modèle modélisé par la fonction staircase() (à utiliser avec le nls() fonction).
• veuillez noter que j'utilise le purrr y glue dans l'exemple ci-dessous.

Code

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

Toute aide est la bienvenue.

Answer 1

Essayez plutôt le 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

Answer 2

Je suppose que l'on vous donne un vecteur d'observations de longueur len comme celles tracées dans votre exemple, et vous souhaitez identifier k sauts et k sauter des tailles. (Ou peut-être que je vous ai mal compris ; mais vous n'avez pas vraiment dit ce que vous voulez obtenir). Ci-dessous, je vais esquisser une solution en utilisant la recherche locale. Je commence avec les données de votre exemple :

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)

Une solution est une liste de postes y tailles des sauts. Notez que j'utilise des vecteurs pour stocker ces données, car il devient fastidieux de définir des variables lorsque vous avez 20 sauts, par exemple.

Un exemple de solution (aléatoire) :

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

Nous avons besoin d'une fonction objective pour calculer la qualité de la solution. Je définis également une fonction d'aide simple stairs qui est utilisé par la fonction objective. La fonction objectif abs_diff calcule la différence absolue moyenne entre la série ajustée (telle que définie par la solution) et 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)
}

Vient maintenant le composant clé de la recherche locale : la fonction de voisinage qui est utilisée pour faire évoluer la solution. La fonction de voisinage prend une solution et la modifie légèrement. Ici, elle va soit choisir une position ou un taille et le modifier légèrement.

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

Un exemple d'appel : ici, la taille du premier saut de la nouvelle solution est modifiée.

## > 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

Je reste pour lancer la recherche locale.

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

Nous pouvons tracer la solution : la ligne ajustée est représentée en bleu.

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