Pourquoi nls () me donne-t-il des erreurs de «matrice de gradient singulière aux estimations initiales des paramètres»?

21

J'ai quelques données de base sur les réductions d'émissions et le coût par voiture:

q24 <- read.table(text = "reductions  cost.per.car
    50  45
    55  55
    60  62
    65  70
    70  80
    75  90
    80  100
    85  200
    90  375
    95  600
    ",header = TRUE, sep = "")

Je sais que c'est une fonction exponentielle, donc je m'attends à pouvoir trouver un modèle qui correspond à:

    model <- nls(cost.per.car ~ a * exp(b * reductions) + c, 
         data = q24, 
         start = list(a=1, b=1, c=0))

mais je reçois une erreur:

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

J'ai lu une tonne de questions sur l'erreur que je vois et je comprends que le problème est probablement que j'ai besoin de startvaleurs meilleures / différentes (cela initial parameter estimatesfait un peu plus de sens) mais je ne suis pas sûr, étant donné le les données dont je dispose, comment je procéderais pour estimer de meilleurs paramètres.

Amanda
la source
Je suggère de commencer votre déchiffrement en recherchant sur notre site le message d'erreur .
whuber
3
En fait, je l'ai fait et ma recherche de l'erreur complète a révélé une question à moitié cuite avec trois points de données et aucune réponse. Mais votre recherche plus spécifique donne des résultats. Peut-être parce que vous avez plus d'expérience ici et que vous savez quels termes se distinguent comme pertinents.
Amanda
Une chose que j'ai trouvée au sujet des erreurs logicielles est qu'une recherche du message d'erreur spécifique (généralement entre guillemets) est le moyen le plus sûr de savoir s'il a été discuté auparavant. (Cela vaut pour Internet, pas seulement sur les sites SE.) Comme le dit notre message "en attente", si vos recherches supplémentaires ne résolvent pas votre problème, veuillez revenir et nous repousser un peu: cette question est à l'intersection des statistiques et de l'informatique et pourrait exposer ici quelques questions d'un grand intérêt.
whuber
1
L'ajustement pour vos valeurs de départ est très loin des données; comparer exp(50)et exp(95)aux valeurs y à x = 50 et x = 95. Si vous définissez c=0et prenez un log de y (établissant une relation linéaire), vous pouvez utiliser la régression pour obtenir des estimations initiales pour log ( ) et b qui suffiront pour vos données (ou si vous ajustez une ligne à travers l'origine, vous pouvez quitter a à 1 et utilisez simplement l'estimation pour b ; cela suffit également pour vos données). Si b est bien en dehors d'un intervalle assez étroit autour de ces deux valeurs, vous rencontrerez des problèmes. [Alternativement, essayez un autre algorithme]ababb
Glen_b -Reinstate Monica
1
Merci @Glen_b. J'espérais pouvoir utiliser R au lieu d'une calculatrice graphique pour travailler à travers un manuel d'introduction aux statistiques (et sauter le cours lui-même), donc je commence avec seulement les informations statistiques les plus élémentaires, mais beaucoup d'expérience avec d'autres découpages et découpages en R .
Amanda

Réponses:

38

Trouver automatiquement de bonnes valeurs de départ pour un modèle non linéaire est un art. (Il est relativement facile pour des ensembles de données uniques lorsque vous pouvez simplement tracer les données et faire de bonnes suppositions visuellement.) Une approche consiste à linéariser le modèle et à utiliser les estimations des moindres carrés.

Dans ce cas, le modèle a la forme

E(Y)=aexp(bx)+c

pour les paramètres inconnus . La présence de l'exponentielle nous encourage à utiliser des logarithmes - mais l'ajout de c rend cela difficile. Remarquez, cependant, que si un est alors positif c sera inférieure à la plus petite valeur attendue de Y --et pourrait donc être un peu moins que la plus petite observée valeur de Y . (Si a peut être négatif, vous devrez également considérer une valeur dea,b,ccacYYa légèrement supérieure à la plus grande valeur observée de Y. )cY

Prenons donc soin de en utilisant comme estimation initiale c 0 quelque chose comme la moitié du minimum des observationscc0 . Le modèle peut maintenant être réécrit sans ce terme additif épineux commeyi

E(Y)c0aexp(bx).

Que nous pouvons prendre le journal de:

log(E(Y)c0)log(a)+bx.

Il s'agit d'une approximation linéaire du modèle. Les deux log(a) et peuvent être estimés avec moindres carrés.b

Voici le code révisé:

c.0 <- min(q24$cost.per.car) * 0.5
model.0 <- lm(log(cost.per.car - c.0) ~ reductions, data=q24)
start <- list(a=exp(coef(model.0)[1]), b=coef(model.0)[2], c=c.0)
model <- nls(cost.per.car ~ a * exp(b * reductions) + c, data = q24, start = start)

Sa sortie (pour les données d'exemple) est

Nonlinear regression model
  model: cost.per.car ~ a * exp(b * reductions) + c
   data: q24
        a         b         c 
 0.003289  0.126805 48.487386 
 residual sum-of-squares: 2243

Number of iterations to convergence: 38 
Achieved convergence tolerance: 1.374e-06

La convergence semble bonne. Trouvons-le:

plot(q24)
p <- coef(model)
curve(p["a"] * exp(p["b"] * x) + p["c"], lwd=2, col="Red", add=TRUE)

Figure

Cela a bien fonctionné!

Lors de l'automatisation, vous pouvez effectuer des analyses rapides des résidus, comme comparer leurs extrêmes à la dispersion dans les données ( ). Vous pourriez également avoir besoin d'un code analogue pour faire face à la possibilité d' un < 0ya<0 ; Je laisse cela comme un exercice.


Une autre méthode pour estimer les valeurs initiales repose sur la compréhension de leur signification, qui peut être basée sur l'expérience, la théorie physique, etc. Un exemple étendu d'un ajustement non linéaire (modérément difficile) dont les valeurs initiales peuvent être déterminées de cette manière est décrit dans ma réponse. à /stats//a/15769 .

Analyse visuelle d'un nuage de points (pour déterminer les estimations initiales des paramètres) est décrite et illustrée sur /stats//a/32832 .

Dans certaines circonstances, une séquence d'ajustements non linéaires est effectuée où vous pouvez vous attendre à ce que les solutions changent lentement. Dans ce cas, il est souvent pratique (et rapide) d' utiliser les solutions précédentes comme estimations initiales pour les suivantes . Je me souviens avoir utilisé cette technique (sans commentaire) sur /stats//a/63169 .

Whuber
la source
2

Cette bibliothèque a pu résoudre mon problème avec les nls singular gradient: http://www.r-bloggers.com/a-better-nls/ Un exemple:

library(minpack.lm)
nlsLM(function, start=list(variable=2,variable2=12))
joshring
la source
Cette fonction semble être appelée nls.lmmaintenant.
Matt
-1

Donc ... je pense que j'ai mal interprété cela comme une fonction exponentielle. Tout ce dont j'avais besoin étaitpoly()

model <- lm(cost.per.car ~ poly(reductions, 3), data=q24)
new.data <- data.frame(reductions = c(91,92,93,94))
predict(model, new.data)

plot(q24)
lines(q24$reductions, predict(model, list(reductions = q24$reductions)))

Ou en utilisant lattice:

xyplot(cost.per.car ~ reductions, data = q24,
       panel = function(x, y) {
         panel.xyplot(x, y)
         panel.lines(x, predict(model,list(reductions = x) ))
       }, 
       xlab = "Reductions", 
       ylab = "Cost per car")
Amanda
la source
2
Cela ne répond pas à la question que vous avez posée - cela la change en quelque chose de différent (et plutôt moins intéressant, à mon humble avis).
whuber
6
Bien que cela puisse résoudre le problème de l'ajustement d'une fonction pour représenter les données, vos réponses acceptées ne sont pas attendues pour votre question. M. @whuber vous a fourni une excellente explication et mérite la réponse acceptée.
Lourenco