<div dir="ltr"><div>Pessoal, boa noite!</div><div><br></div><div>Tenho um conjunto de dados que provavelmente segue uma distribuição Weibull. </div><div><br></div><div>Os dados são os seguintes:</div><div><br></div><div>y<-c(295000,</div><div>869000<span class="" style="white-space:pre">       </span>,</div><div>869900<span class="" style="white-space:pre">    </span>,</div><div>1573335<span class="" style="white-space:pre">   </span>,</div><div>151400<span class="" style="white-space:pre">    </span>,</div><div>152000<span class="" style="white-space:pre">    </span>,</div><div>183700<span class="" style="white-space:pre">    </span>,</div><div>218000<span class="" style="white-space:pre">    </span>,</div><div>30200<span class="" style="white-space:pre">     </span>,</div><div>45100<span class="" style="white-space:pre">     </span>,</div><div>46900<span class="" style="white-space:pre">     </span>,</div><div>47300<span class="" style="white-space:pre">     </span></div><div>)</div><div><br></div><div><br></div><div>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:</div><div><br></div><div>###### Newton Paphson Method ##########</div><div><br></div><div>n=length(y)</div><div>x=vector()</div><div><br></div><div>##### Função que quero obter a raiz (estimador de lambda)</div><div>f <- quote((sum((y^delta)*log(y))/sum(y^delta))-(1/delta)-(1/n)*sum(log(y)))</div><div>fx0 <- function(delta){ eval(f) }</div><div><br></div><div>##### Primeira derivada da função (estimador de lambda)</div><div>f1 <-quote((1/delta^2)-((sum((log(y))*(y^delta)))^2/(sum(y^delta))^2)+(sum((log(y))^2*(y^delta))/sum(y^delta)))</div><div>fx1 <- function(delta){ eval(f1) }</div><div><br></div><div>### Obtendo a estimativa para delta</div><div><span class="" style="white-space:pre">  </span>i <- 1       # valor inicial para o passo</div><div><span class="" style="white-space:pre">       </span>dif <- 1   # valor inical para a diferença entre duas avaliações</div><div><span class="" style="white-space:pre">        </span>tol <- 10^-6 # diferência mínima entre duas avaliações (tolerância)</div><div><span class="" style="white-space:pre">     </span>passo <- 1000   # número máximo de passos</div><div><span class="" style="white-space:pre">       </span>x <- 0.3      # valor inicial para a raiz<span class="" style="white-space:pre">      </span>             <span class="" style="white-space:pre">     </span></div><div><span class="" style="white-space:pre">           </span>while(i<passo & dif>tol){</div><div>  <span class="" style="white-space:pre">              </span>x[i+1] <- x[i]-(fx0(x[i])/fx1(x[i]))</div><div>  <span class="" style="white-space:pre">          </span>dif <- abs(x[i+1]-x[i])</div><div><span class="" style="white-space:pre">         </span>d<-x[i] </div><div><span class="" style="white-space:pre">                </span>i<-i+1</div><div>}</div><div><br></div><div>### Obtendo a estimativa para alpha</div><div><br></div><div>alpha=((sum(y^d))/n)^(1/d)</div><div><br></div><div>### d é o parâmetro de forma e alpha é o parâmetro de escala</div><div><br></div><div>d</div><div>alpha</div><div><br></div><div><br></div><div>########## Teste de aderência </div><div><br></div><div>ypop<-rweibull(1000,d,alpha)</div><div>ks.test(y,ypop)</div><div><br></div><div>#####################################################</div><div><br></div><div><br></div><div><br></div><div>###### Nelder-Mead através da função Optim ##########</div><div><br></div><div><br></div><div>n=length(y)</div><div>vero <- function(par,x){</div><div>k = par[1] #### Shape parameter</div><div>lambda = par[2] #### Scale Parameter</div><div>   saida<-(n*log(k/lambda) + (k-1)*sum(log(x/lambda)) - sum((x/lambda)^k))</div><div>return(-saida)</div><div>}</div><div><br></div><div>saiday<-optim(par=c(1.8,4),fn=vero,</div><div>      method= "Nelder-Mead",x=y</div><div>      )</div><div>k<-saiday[1]$par[1]</div><div>lambda<-saiday[1]$par[2]</div><div><br></div><div>### k é o parâmetro de forma e lambda e o parâmetro de escala</div><div><br></div><div>k</div><div>lambda</div><div><br></div><div>########## Teste de aderência </div><div><br></div><div>ypop<-rweibull(1000,k,lambda)</div><div>ks.test(y,ypop)</div><div><br></div><div><br></div><div>###############################################</div><div><br></div><div><br></div><div><br></div><div>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....</div><div><br></div><div><br></div><div>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???</div><div><br></div><div><br></div><div><br></div><div>Observe que se eu gerar uma weibull, por exemplo:</div><div><br></div><div>Y<-rweibull(12,5,6)</div><div><br></div><div>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)...</div><div><br></div><div>O que pode estar acontecendo?? Vcs tem alguma idéia??? O problema está nos dados que estou analizando?? que tipo de problema seria esse??</div><div><br></div><div><br></div><div>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... </div><div><br></div><div><br></div><div>Agradeço qualquer dica!</div><div><br></div><div>Abraço,</div><div><br></div><div>Romero.</div><div><br></div><br><div class="gmail_signature"><div></div></div>
</div>