[R-br] Ajustamento de dados à curva exponencial

Quinta Novembro 28 16:11:57 BRST 2013

Caro Walmes, boa tarde.

Li o comentário sobre “curva exponencial” no e-mail abaixo endereçado ao
grupo.

Pois bem, estou com um probleminha para resolver, a saber:

Tenho estas duas séries:

x<-c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23)

y<-c(100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,950
,1000,1050,1100,1150,1200)

Utilizei um script apresentado no Livro Conhecendo o R de Luiz A Pertenelli
e Márcio P de Mello pág. 165 e não consegui resolver no R.

Bom! Ajustei a curva y para exponencial Y no EXCEL conforme abaixo:

Y<-c(187,206,227,250,275,303,334,368,405,446,491,541,596,657,723,797,877,966
,1064,1172,1291,1422,1566)

Por favor, como ajusto à exponencial a série y no R? Como chegar ao
resultado Y no R?

#Script utilizado que não funcionou:

x<-c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23)

y<-c(100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,950
,1000,1050,1100,1150,1200)

dados<-data.frame(x,y)

plot(dados)

funcao<-y~a*exp(b*x)

exponencial<-nls(funcao, dados, start=c(a=1,b=1))

exponencial

=====Erro no R ==================

10,11,12,13,14,15,16,17,18,19,20,21,22,23)

>
y<-c(100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,950
,1000,1050,1100,1150,1200)

> dados<-data.frame(x,y)

> plot(dados)

> funcao<-y~a*exp(b*x)

> exponencial<-nls(funcao, dados, start=c(a=1,b=1))

Erro em nls(funcao, dados, start = c(a = 1, b = 1)) : gradiente singular

> exponencial

Erro: objeto 'exponencial' não encontrado

> 

Onde estou errando?

Abs.

  _____  

 <http://www.posicional.com/> assinatura_email_sites

Cícero C Nunes

Diretor Técnico Operacional

Divisão Remuneração

 <mailto:c2n em posicional.com> c2n em posicional.com

+55(16)3397.0226

+55(11)3280.0226

+55(11)9.9978.0026

SKYPE: consultormp

Rua Padre Duarte, 1360 – salas 6/7 – Centro – Araraquara, SP – CEP.:
14.801-310

“Antes de imprimir pense em sua responsabilidade com o meio ambiente – Quem
ama cuida!”

De: R-br [mailto:r-br-bounces em listas.c3sl.ufpr.br] Em nome de walmes .
Enviada em: sexta-feira, 8 de novembro de 2013 18:33
Para: r-br em listas.c3sl.ufpr.br
Assunto: Re: [R-br] Contrastes em modelos de regressão de poisson
inflacionados por zeros

Você não deve esperar ver uma curva exponencial para a média das contagens
em função do tempo e uma sigmoide para a proporção de zeros não Poisson para
os zeros se não simulou os dados com efeito de tempo e não declarou efeito
de tempo na porção binomial. Para que isso aconteça a simulação tem que
estar condizente, tem que ter os efeitos para que as curvas não sejam as
"curvas nulas".

#------------------------------------------------------------------
# Definições da sessão.

rm(list=ls())
require(pscl)
require(VGAM)
require(multcomp)
require(lattice)
require(latticeExtra)

#------------------------------------------------------------------
# Dados artificiais.

da <- expand.grid(tempo=rep(1:10), trat=gl(3,10))

## simula porção Poisson
X <- model.matrix(~trat+tempo, da)
colnames(X)
beta.pois <- c(-1,0.25,0.5,0.25)
eta.pois <- X%*%beta.pois
da$y.pois <- rpois(nrow(X), lambda=exp(eta.pois))
xyplot(y.pois~tempo|trat, da)

## simula porção binomial
X <- model.matrix(~tempo, da)
beta.bern <- c(0,0.3)
eta.bern <- X%*%beta.bern
da$y.bern <- rbinom(nrow(X), size=1, prob=1/(1+exp(-eta.bern)))

xyplot(y.bern~tempo|trat, da)

## monta a contagem com inflação de zeros
da$y <- with(da, ifelse(y.bern==0, 0, y.pois))
xyplot(y~tempo|trat, da)

#------------------------------------------------------------------
# Modelo completo.

compl.mod <- zeroinfl(y~trat+tempo|trat+tempo, data=da)
coef(compl.mod)

#------------------------------------------------------------------
# Predição do modelo considerando as duas porções.

pred <- expand.grid(tempo=rep(1:10), trat=gl(3,1))
X <- model.matrix(~trat+tempo, pred)
i <- grep("^count\\_", names(coef(compl.mod)))
eta <- X%*%coef(compl.mod)[i]
pred$m.pois <- exp(eta)
i <- grep("^zero\\_", names(coef(compl.mod)))
eta <- X%*%coef(compl.mod)[i]
pred$p.zero <- exp(eta)/(1+exp(eta))

xyplot(y~tempo|trat, data=da, jitter.x=TRUE)+
    as.layer(xyplot(m.pois~tempo|trat, data=pred, type="l"))+
    as.layer(xyplot(p.zero~tempo|trat, data=pred,
                    type="l", lty=2, lwd=2))+
    layer(panel.abline(h=1, lty=2))

xyplot(p.zero~tempo|trat, data=pred,
                    type="l", lty=2, lwd=2)

#------------------------------------------------------------------

À disposição.
Walmes.

==========================================================================
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
skype: walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes <http://www.leg.ufpr.br/%7Ewalmes> 
linux user number: 531218
==========================================================================

-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20131128/7664b9dd/attachment.html>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: image001.gif
Tipo: image/gif
Tamanho: 12852 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20131128/7664b9dd/attachment.gif>