Estimativas de parâmetros de uma Weibull utilizando duas metodologias: Newton-Raphson "vs" Nelder Mead (estimativas muito diferentes) o que pode estar acontecendo??
Pessoal, boa noite! Tenho um conjunto de dados que provavelmente segue uma distribuição Weibull. Os dados são os seguintes: y<-c(295000, 869000 , 869900 , 1573335 , 151400 , 152000 , 183700 , 218000 , 30200 , 45100 , 46900 , 47300 ) Efetuei a estimação dos parâmetros de duas formas, uma utilizando as expressões analíticas onde tive que usar o método de Newton-Raphson... E outra utilizando a função optim em que selecionei o método de Nelder Mead... Os dois algoritmos seguem abaixo: ###### Newton Paphson Method ########## n=length(y) x=vector() ##### Função que quero obter a raiz (estimador de lambda) f <- quote((sum((y^delta)*log(y))/sum(y^delta))-(1/delta)-(1/n)*sum(log(y))) fx0 <- function(delta){ eval(f) } ##### Primeira derivada da função (estimador de lambda) f1 <-quote((1/delta^2)-((sum((log(y))*(y^delta)))^2/(sum(y^delta))^2)+(sum((log(y))^2*(y^delta))/sum(y^delta))) fx1 <- function(delta){ eval(f1) } ### Obtendo a estimativa para delta i <- 1 # valor inicial para o passo dif <- 1 # valor inical para a diferença entre duas avaliações tol <- 10^-6 # diferência mínima entre duas avaliações (tolerância) passo <- 1000 # número máximo de passos x <- 0.3 # valor inicial para a raiz while(i<passo & dif>tol){ x[i+1] <- x[i]-(fx0(x[i])/fx1(x[i])) dif <- abs(x[i+1]-x[i]) d<-x[i] i<-i+1 } ### Obtendo a estimativa para alpha alpha=((sum(y^d))/n)^(1/d) ### d é o parâmetro de forma e alpha é o parâmetro de escala d alpha ########## Teste de aderência ypop<-rweibull(1000,d,alpha) ks.test(y,ypop) ##################################################### ###### Nelder-Mead através da função Optim ########## n=length(y) vero <- function(par,x){ k = par[1] #### Shape parameter lambda = par[2] #### Scale Parameter saida<-(n*log(k/lambda) + (k-1)*sum(log(x/lambda)) - sum((x/lambda)^k)) return(-saida) } saiday<-optim(par=c(1.8,4),fn=vero, method= "Nelder-Mead",x=y ) k<-saiday[1]$par[1] lambda<-saiday[1]$par[2] ### k é o parâmetro de forma e lambda e o parâmetro de escala k lambda ########## Teste de aderência ypop<-rweibull(1000,k,lambda) ks.test(y,ypop) ############################################### O que acontece é que as estimativas através dos dois métodos estão dando muito diferentes e através de um KS.test concluo que as estimativas através de Newton-Raphson são muito boas, pois o p-valor fica realmente bem alto; já as estimativas utilizando Nelder-Mead ficam péssimas e o ajuste, claro, fica péssimo com p-valor muito pequeno.... Gostaria de saber se existe alguma explicação para que isso esteja ocorrendo... O método de Nelder-Mead não é totalmente indicado?? O que estaria provocando essa disparidade nas estimativas... Sei que são dois métodos de otimização diferentes, mas acredito que os dois eram pra fornecer vaores bem próximos, concordam??? Observe que se eu gerar uma weibull, por exemplo: Y<-rweibull(12,5,6) As estimativas serão bem próximas e o ks.test me dará um p-valor bem condizente com a realidade (isso para os dois métodos de otimização)... O que pode estar acontecendo?? Vcs tem alguma idéia??? O problema está nos dados que estou analizando?? que tipo de problema seria esse?? Estou confuso pois quando gero dados no R os dois métodos apresentam boas estimtivas e p-valor do ks.test é bem condizente... Entretanto se uso os dados que tenho, os parâmetros só são bem estimados pelo método de Newton-Raphson o outro método fica péssimo... Agradeço qualquer dica! Abraço, Romero.
Provavelmente seu Nelder-Mead não convergiu ou convergiu para algum minimo local. Vc tentou usar valores iniciais diferentes? Pra sair da dúvida pegue as estimativas obtivas pelo newton-raphson e coloque como inicial do Nelder-Mead e veja o que acontece. Se o Nelder-Mead não convergir para o mesmo valor vc pode ter algum erro na função de log-verossimilhança. E tbm vc pode tentar os outros algorithmos disponiveis na optim e por fim use o pacote bbmle para fazer os perfis de verossimilhança além da informação de como são os perfis isso ainda funciona como um check se seu algoritmo realmente convergiu pro lugar certo. Em 7 de maio de 2015 00:20, Romero Luiz M. Sales Filho < romero.sfilho@gmail.com> escreveu:
Pessoal, boa noite!
Tenho um conjunto de dados que provavelmente segue uma distribuição Weibull.
Os dados são os seguintes:
y<-c(295000, 869000 , 869900 , 1573335 , 151400 , 152000 , 183700 , 218000 , 30200 , 45100 , 46900 , 47300 )
Efetuei a estimação dos parâmetros de duas formas, uma utilizando as expressões analíticas onde tive que usar o método de Newton-Raphson... E outra utilizando a função optim em que selecionei o método de Nelder Mead... Os dois algoritmos seguem abaixo:
###### Newton Paphson Method ##########
n=length(y) x=vector()
##### Função que quero obter a raiz (estimador de lambda) f <- quote((sum((y^delta)*log(y))/sum(y^delta))-(1/delta)-(1/n)*sum(log(y))) fx0 <- function(delta){ eval(f) }
##### Primeira derivada da função (estimador de lambda) f1 <-quote((1/delta^2)-((sum((log(y))*(y^delta)))^2/(sum(y^delta))^2)+(sum((log(y))^2*(y^delta))/sum(y^delta))) fx1 <- function(delta){ eval(f1) }
### Obtendo a estimativa para delta i <- 1 # valor inicial para o passo dif <- 1 # valor inical para a diferença entre duas avaliações tol <- 10^-6 # diferência mínima entre duas avaliações (tolerância) passo <- 1000 # número máximo de passos x <- 0.3 # valor inicial para a raiz while(i<passo & dif>tol){ x[i+1] <- x[i]-(fx0(x[i])/fx1(x[i])) dif <- abs(x[i+1]-x[i]) d<-x[i] i<-i+1 }
### Obtendo a estimativa para alpha
alpha=((sum(y^d))/n)^(1/d)
### d é o parâmetro de forma e alpha é o parâmetro de escala
d alpha
########## Teste de aderência
ypop<-rweibull(1000,d,alpha) ks.test(y,ypop)
#####################################################
###### Nelder-Mead através da função Optim ##########
n=length(y) vero <- function(par,x){ k = par[1] #### Shape parameter lambda = par[2] #### Scale Parameter saida<-(n*log(k/lambda) + (k-1)*sum(log(x/lambda)) - sum((x/lambda)^k)) return(-saida) }
saiday<-optim(par=c(1.8,4),fn=vero, method= "Nelder-Mead",x=y ) k<-saiday[1]$par[1] lambda<-saiday[1]$par[2]
### k é o parâmetro de forma e lambda e o parâmetro de escala
k lambda
########## Teste de aderência
ypop<-rweibull(1000,k,lambda) ks.test(y,ypop)
###############################################
O que acontece é que as estimativas através dos dois métodos estão dando muito diferentes e através de um KS.test concluo que as estimativas através de Newton-Raphson são muito boas, pois o p-valor fica realmente bem alto; já as estimativas utilizando Nelder-Mead ficam péssimas e o ajuste, claro, fica péssimo com p-valor muito pequeno....
Gostaria de saber se existe alguma explicação para que isso esteja ocorrendo... O método de Nelder-Mead não é totalmente indicado?? O que estaria provocando essa disparidade nas estimativas... Sei que são dois métodos de otimização diferentes, mas acredito que os dois eram pra fornecer vaores bem próximos, concordam???
Observe que se eu gerar uma weibull, por exemplo:
Y<-rweibull(12,5,6)
As estimativas serão bem próximas e o ks.test me dará um p-valor bem condizente com a realidade (isso para os dois métodos de otimização)...
O que pode estar acontecendo?? Vcs tem alguma idéia??? O problema está nos dados que estou analizando?? que tipo de problema seria esse??
Estou confuso pois quando gero dados no R os dois métodos apresentam boas estimtivas e p-valor do ks.test é bem condizente... Entretanto se uso os dados que tenho, os parâmetros só são bem estimados pelo método de Newton-Raphson o outro método fica péssimo...
Agradeço qualquer dica!
Abraço,
Romero.
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- Wagner Hugo Bonat ---------------------------------------------------------------------------------------------- Department of Mathematics and Computer Science (IMADA) University of Southern Denmark (SDU) and Laboratório de Estatística e Geoinformação (LEG) Universidade Federal do Paraná (UFPR)
que tal um bom chute inicial? hist(y, freq=FALSE) plot(function(x) dweibull(x, 1, 1e6), 0, 2e6, add=T, col=1) plot(function(x) dweibull(x, 1.2, 1e6), 0, 2e6, add=T, col=2) essa aplicacao do ks.test() parece que nao faz muito sentido aqui...
seria bom mudar a escala dos dados para evitar problemas numericos ### estima para y/1000000 o <- optim(c(1.2, 1), function(p) -sum(dweibull(y/1e6, p[1], p[2], log=TRUE)), method='BFGS', hessian=TRUE) ### estima para 1000 reamostras (bootstrap) oo <- t(replicate(1000, optim(o$par, function(p) -sum(dweibull(sample(y, replace=TRUE)/1e6, p[1], p[2], log=TRUE)), method='BFGS')$par)) (tb <- rbind(sd.assintotico=sqrt(diag(solve(o$hess))), sd.bootstrap=apply(oo,2,sd))) ### visualiza estimativas bootstrap e adiciona SD assintotico par(mfrow=c(2,2), mar=c(3,3,1,1), mgp=c(2,1,0)) hist(y/1e6, main='dados') plot(function(x) dweibull(x, o$par[1], o$par[2]), 0, 2, add=T, col=2) hist(oo[,1], main='shape') abline(v=o$par[1]+c(-2,0,2)*tb[1,1], col=2, lty=c(2,1,2)) hist(oo[,2], main='scale') abline(v=o$par[2]+c(-2,0,2)*tb[1,2], col=2, lty=c(2,1,2)) Elias
participantes (3)
-
Elias T. Krainski -
Romero Luiz M. Sales Filho -
Wagner Bonat