Como minimizar a soma residual dos quadrados de um ajuste exponencial?

14

Eu tenho os seguintes dados e gostaria de ajustar um modelo de crescimento exponencial negativo a ele:

Days <- c( 1,5,12,16,22,27,36,43)
Emissions <- c( 936.76, 1458.68, 1787.23, 1840.04, 1928.97, 1963.63, 1965.37, 1985.71)
plot(Days, Emissions)
fit <- nls(Emissions ~ a* (1-exp(-b*Days)), start = list(a = 2000, b = 0.55))
curve((y = 1882 * (1 - exp(-0.5108*x))), from = 0, to =45, add = T, col = "green", lwd = 4)

O código está funcionando e uma linha de ajuste é plotada. No entanto, o ajuste não é visualmente ideal e a soma residual dos quadrados parece ser bastante grande (147073).

Como podemos melhorar nosso ajuste? Os dados permitem um melhor ajuste?

Não conseguimos encontrar uma solução para esse desafio na rede. Qualquer ajuda direta ou ligação a outros sites / postagens é muito apreciada.

Strohmi
fonte
1
Nesse caso, se você considerar um modelo de regressão , onde , obterá estimadores semelhantes. Ao traçar as regiões de confiança, pode-se observar como esses valores estão contidos nas regiões de confiança. Você não pode esperar um ajuste perfeito, a menos que interpole os pontos ou use um modelo não-linear mais flexível. ϵ iN ( 0 , σ )Emissionsi=f(Daysi,a,b)+ϵiϵiN(0,σ)
Mudei o título porque "modelo exponencial negativo" significa algo diferente do descrito na pergunta.
whuber
Agradecemos por esclarecer a pergunta (@whuber) e por sua resposta (@Procrastinator). Como posso calcular e plotar as regiões de confiança. E qual seria um modelo não-linear mais flexível?
Strohmi
4
Você precisa de um parâmetro adicional. Veja o que acontece com fit <- nls(Emissions ~ a* (1- u*exp(-b*Days)), start = list(a = 2000, b = 0.1, u=.5)); beta <- coefficients(fit); curve((y = beta["a"] * (1 - beta["u"] * exp(-beta["b"]*x))), add = T).
whuber
1
@ whuber - talvez você deva postar isso como resposta?
23412 jbowman

Respostas:

16

Uma lei exponencial (negativa) assume a forma . Ao permitir alterações de unidades dos x e y valores, no entanto, dizer para y = α y ' + β e x = γ x ' + δ , então a lei será expressa comoy=exp(x)xyy=αy+βx=γx+δ

αy+β=y=exp(x)=exp(γxδ),

que algebricamente é equivalente a

y=1αexp(γxδ)β=a(1uexp(bx))

utilizando-se três parâmetros de , u = 1 / ( β exp ( δ ) ) , e b = γ . Podemos reconhecer a como um parâmetro de escala para y , b como um parâmetro de escala para x e u como derivado de um parâmetro de localização para x .a=β/αu=1/(βexp(δ))b=γaybxux

Como regra geral, esses parâmetros podem ser identificados rapidamente na plotagem :

  • O parâmetro é o valor da assíntota horizontal, um pouco menos que 2000 .a2000

  • O parâmetro é a quantidade relativa em que a curva aumenta da origem à sua assíntota horizontal. Aqui, o aumento é, portanto, um pouco menos que 2000 - 937 ; relativamente, isso é cerca de 0,55 da assíntota.u20009370.55

  • Como , quando x é igual a três vezes o valor de 1 / b, a curva deve ter subido para cerca de 1 - 0,05 ou 95 % de seu total. 95 % do aumento de 937 para quase 2000 nos coloca por volta de 1950 ; a varredura no gráfico indica que isso levou de 20 a 25 dias. Chamado de Let It 24 de simplicidade, de onde b 3 / 24exp(3)0.05x1/b10.0595%95%93720001950202524 . (Essemétodo de 95 % para observar uma escala exponencial é padrão em alguns campos que usam muito plotagens exponenciais.)b3/24=0.12595%

Vamos ver como é isso:

plot(Days, Emissions)
curve((y = 2000 * (1 - 0.56 * exp(-0.125*x))), add = T)

Ajuste do globo ocular

Nada mal para um começo! (Mesmo apesar de digitar 0.56no lugar de 0.55, que era uma aproximação grosseira de qualquer maneira.) Podemos polir isso com nls:

fit <- nls(Emissions ~ a * (1- u * exp(-b*Days)), start=list(a=2000, b=1/8, u=0.55))
beta <- coefficients(fit)
plot(Days, Emissions)
curve((y = beta["a"] * (1 - beta["u"] * exp(-beta["b"]*x))), add = T, col="Green", lwd=2)

Ajuste NLS

A saída de nlscontém informações abrangentes sobre incerteza de parâmetro. Por exemplo , um simples summaryfornece erros padrão de estimativas:

> summary(fit)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
a 1.969e+03  1.317e+01  149.51 2.54e-10 ***
b 1.603e-01  1.022e-02   15.69 1.91e-05 ***
u 6.091e-01  1.613e-02   37.75 2.46e-07 ***

Podemos ler e trabalhar com toda a matriz de covariância das estimativas, o que é útil para estimar intervalos de confiança simultâneos (pelo menos para grandes conjuntos de dados):

> vcov(fit)
             a             b             u
a 173.38613624 -8.720531e-02 -2.602935e-02
b  -0.08720531  1.044004e-04  9.442374e-05
u  -0.02602935  9.442374e-05  2.603217e-04

nls suporta gráficos de perfil para os parâmetros, fornecendo informações mais detalhadas sobre suas incertezas:

> plot(profile(fit))

a

Gráfico de perfil

219451995

whuber
fonte
res <- residuals(fit); res %*% resu2724147073
Tudo bem e bom whuber. Mas talvez o OP tivesse algum motivo para escolher o modelo exponencial (ou talvez seja apenas porque é bem conhecido). Acho que primeiro os resíduos devem ser analisados ​​para o modelo exponencial. Plote-os contra possíveis covariáveis ​​para ver se há estrutura lá e não apenas um grande ruído aleatório. Antes de pular para modelos mais sofisticados, tente ver se um modelo mais sofisticado poderia ajudar.
22912 Michael Michael Chernick
3
x
2
Eu não estava criticando sua resposta! Não vi parcelas residuais. Tudo o que eu estava sugerindo é que gráficos de resíduos versus covariáveis ​​em potencial deveriam ser o primeiro passo para encontrar um modelo melhor. Se eu pensasse que tinha uma resposta para colocar lá em cima, teria dado uma resposta em vez de ter levantado meu argumento como uma constante. Eu pensei que você tivesse dado uma ótima resposta e eu estava entre os que lhe deram um +1.
22812 Michael R. Chernick