Por que nls () está me dando erros de “matriz de gradiente singular nas estimativas de parâmetros iniciais”?

21

Eu tenho alguns dados básicos sobre reduções de emissões e custo por carro:

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

Sei que essa é uma função exponencial, portanto, espero encontrar um modelo que se encaixe com:

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

mas estou recebendo um erro:

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

Eu li várias perguntas sobre o erro que estou vendo e estou percebendo que o problema provavelmente é que preciso de startvalores melhores / diferentes (o que initial parameter estimatesfaz um pouco mais de sentido), mas não tenho certeza, considerando o dados que tenho, como calcularia melhores parâmetros.

Amanda
fonte
Sugiro que você inicie sua decifração pesquisando em nosso site a mensagem de erro .
whuber
3
Na verdade, fiz isso e minha busca pelo erro completo resultou em uma pergunta incompleta com três pontos de dados e nenhuma resposta. Mas sua pesquisa mais específica obtém alguns resultados. Possivelmente porque você tem mais experiência aqui e sabe quais termos se destacam como relevantes.
Amanda
Uma coisa que eu descobri sobre erros de software é que uma busca pela mensagem de erro específica (geralmente entre aspas) é a maneira mais segura de descobrir se ela já foi discutida anteriormente. (Isso é válido para toda a Internet, e não apenas para sites de SE). Como diz a mensagem "em espera", se sua pesquisa adicional não resolver o problema, volte e pressione-nos um pouco: esta pergunta está em a interseção de estatística e computação e pode expor alguns problemas de grande interesse aqui.
whuber
1
A adequação aos seus valores iniciais está muito longe dos dados; compare exp(50)e exp(95)com os valores de y em x = 50 ex = 95. Se você definir c=0e registrar o log de y (estabelecendo um relacionamento linear), poderá usar a regressão para obter estimativas iniciais do log ( ) que serão suficientes para seus dados (ou se você ajustar uma linha na origem, poderá sair em 1 e use apenas a estimativa para ; isso também é suficiente para seus dados). Se estiver muito fora de um intervalo bastante estreito em torno desses dois valores, você terá alguns problemas. [Como alternativa, tente um algoritmo diferente]b a b bababb
Glen_b -Reinstar Monica
1
Obrigado @Glen_b. Eu esperava poder usar R no lugar de uma calculadora gráfica para trabalhar com um livro de introdução de estatísticas (e ultrapassar o curso em si), por isso estou começando apenas com o menor insight estatístico, mas com muita experiência em fatiar e cortar dados em R #
219 Amanda Amanda

Respostas:

38

Encontrar automaticamente bons valores iniciais para um modelo não linear é uma arte. (É relativamente fácil para conjuntos de dados pontuais quando você pode apenas plotar os dados e fazer algumas boas suposições visualmente.) Uma abordagem é linearizar o modelo e usar estimativas de mínimos quadrados.

Nesse caso, o modelo tem a forma

E(Y)=aexp(bx)+c

para parâmetros desconhecidos . A presença do exponencial nos encoraja a usar logaritmos - mas a adição de torna difícil fazer isso. Observe, porém, que, se é positivo, então será menor do que o menor valor esperado de --e, portanto, pode ser um pouco menos do que o menor observado valor de . (Se puder ser negativo, você também terá que considerar um valor de um pouco maior que o maior valor observado de ).c a c Y Y aa,b,ccacYYacY

Vamos, então, cuidar de usando como estimativa inicial algo como metade do mínimo das observações . O modelo agora pode ser reescrito sem esse termo aditivo espinhoso comocc0yi

E(Y)c0aexp(bx).

Que podemos pegar o log de:

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

Essa é uma aproximação linear ao modelo. Ambos e podem ser estimados por mínimos quadrados.log(a)b

Aqui está o código revisado:

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)

Sua saída (para os dados de exemplo) é

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

A convergência parece boa. Vamos traçar isso:

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

Figura

Funcionou bem!

Ao automatizar isso, você pode executar algumas análises rápidas dos resíduos, como comparar seus extremos com a dispersão nos dados ( ). Você também pode precisar de código análogo para lidar com a possibilidade ; Deixo isso como um exercício.a < 0ya<0


Outro método para estimar valores iniciais baseia-se na compreensão do que eles significam, que podem ser baseados em experiência, teoria física etc. Um exemplo extenso de um ajuste não-linear (moderadamente difícil) cujos valores iniciais podem ser determinados dessa maneira é descrito em minha resposta em /stats//a/15769 .

A análise visual de um gráfico de dispersão (para determinar as estimativas iniciais dos parâmetros) é descrita e ilustrada em /stats//a/32832 .

Em algumas circunstâncias, é feita uma sequência de ajustes não lineares em que você pode esperar que as soluções mudem lentamente. Nesse caso, geralmente é conveniente (e rápido) usar as soluções anteriores como estimativas iniciais para as próximas . Lembro-me de usar essa técnica (sem comentários) em /stats//a/63169 .

whuber
fonte
2

Esta biblioteca conseguiu resolver meu problema com o nls singular gradient: http://www.r-bloggers.com/a-better-nls/ Um exemplo:

library(minpack.lm)
nlsLM(function, start=list(variable=2,variable2=12))
Joshring
fonte
Essa função parece ser chamada nls.lmagora.
Matt
-1

Então ... acho que interpretei errado isso como uma função exponencial. Tudo que eu precisava erapoly()

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, usando 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
fonte
2
Isso não responde à pergunta que você fez - muda para algo diferente (e um pouco menos interessante, IMHO).
whuber
6
Embora possa resolver o problema de ajustar uma função para representar os dados, suas respostas aceitas não são a expectativa da sua pergunta. O Sr. @whuber forneceu uma excelente explicação e merece a resposta aceita.
Lourenco