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 start
valeurs meilleures / différentes (cela initial parameter estimates
fait 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.
exp(50)
etexp(95)
aux valeurs y à x = 50 et x = 95. Si vous définissezc=0
et 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]Réponses:
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
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,c c a c Y Y a légèrement supérieure à la plus grande valeur observée de Y. )c Y
Prenons donc soin de en utilisant comme estimation initiale c 0 quelque chose comme la moitié du minimum des observationsc c0 . Le modèle peut maintenant être réécrit sans ce terme additif épineux commeyi
Que nous pouvons prendre le journal de:
Il s'agit d'une approximation linéaire du modèle. Les deuxlog(a) et peuvent être estimés avec moindres carrés.b
Voici le code révisé:
Sa sortie (pour les données d'exemple) est
La convergence semble bonne. Trouvons-le:
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 < 0y a<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 .
la source
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:la source
nls.lm
maintenant.Donc ... je pense que j'ai mal interprété cela comme une fonction exponentielle. Tout ce dont j'avais besoin était
poly()
Ou en utilisant
lattice
:la source