## ## "Limpando a área de trabalho" ## !!! cuidado, vai apagar todos os objetos !!! ## rm(list=ls()) par.ori <- par(no.readonly=TRUE) ## ## Gerando Dados (simulados) ## distribuição exponencial com média 10 ## f(x) = \theta \exp{- \theta x} ; \theta = 1/10 set.seed(2012) dados <- round(rexp(10, rate=1/2), dig=2) dados summary(dados) ## ## Definindo a função de log-verossimilhança para uma observação ## Lik <- function(par, dado, relativa=TRUE){ res <- dexp(dado, rate=par) if(relativa) return(res/dexp(dado, rate=1/dado)) else return(res) } ## gráfico das verosimilhanças (e relativas) individuais de cada observação xmax <- max(1/dados)*3 curve(Lik(x, dado=min(dados), relativa=F), 0.05, xmax, xlab=expression(lambda), ylab=expression(L(lambda)), lty=3) points(dados, rep(0, length(dados))) for(i in dados) curve(Lik(x, dado=i, relativa=F), 0, xmax, add=T, lty=3) curve(Lik(x, dado=min(dados), relativa=T), 0.05, xmax, xlab=expression(lambda), ylab=expression(LR(lambda)), lty=3) points(dados, rep(0, length(dados))) for(i in dados) curve(Lik(x, dado=i, relativa=T), 0, xmax, add=T, lty=3) ## deviances individuais Dev <- function(par, dado){ res <- dexp(dado, rate=par) return(-2*(log(res) - dexp(dado, rate=1/dado, log=TRUE))) } curve(Dev(x, dado=min(dados)), 0.05, xmax, xlab=expression(lambda), ylab=expression(D(lambda)), lty=3) points(dados, rep(0, length(dados))) for(i in dados) curve(Dev(x, dado=i), 0, xmax, add=T, lty=3) ## ## Definindo a função de log-verossimilhança ## (assumindo observações independentes) ## função de log-verossimilhança a partir da densidade lLik <- function(theta, dados){ sapply(theta, function(x) sum(dexp(dados, rate=x, log=TRUE))) } #par(mfrow=c(1,2)) ## fazendo o gráfico da log-verossimilhança #par <- (0:400)/200 curve(lLik(x, dados=dados), 0.05, 2) #plot(par, lLik(par, dados), type="l", xlab=expression(theta), ylab=expression(l(theta))) ## .. de outra forma curve(lLik(theta=x, dados), 0.05, 2, xlab=expression(theta), ylab=expression(l(theta)), col=2) ## Reescrevendo log-verossimilhança em função de n e média ## l(\theta) = n ( \log(\theta) - \theta \bar{y}) lvExp0 <- function(theta, y) { llik <- length(y)*(log(theta) - theta * mean(y)) return(llik) } ## sobrepondo ao gráfico anterior curve(lvExp0(theta=x, y=dados), 0.05, 2, add=T, lty=2, col=4) ## redefinindo novamente, mas agora sem a necessidade de (re)calcular media e n a cada chamada da função ## l(\theta) = n log(\theta) - \theta n \bar{y} lvExp <- function(theta, amostra) { ## amostra = list(n, media) llik <- with(amostra, n* (log(theta) - theta * media)) return(llik) } am <- list(n=length(dados), media=mean(dados)) ## ou, sendo mais cuidadoso com possíveis NA's ... am <- list(n=sum(!is.na(dados)), media=mean(dados, na.rm=TRUE)) curve(lvExp(x, amostra=am), 0.05, 2, col=2) ## "zoom" na função entre 0,10 e 0,40 ## para mostrar o comportamento ao redor do máximo curve(lvExp(x, amostra=am), 0.3, 1) ## Comentários ## - última forma evita cálculos desnecessários ## - verossimilhança assimétrica, mas quase simétrica ao redor do máximo ## ## EMV ## ## Funções de interesse: ## \theta^{r+1} = \theta^{r} - U(\theta) / U'(\theta) ## U(\theta) = l'(\theta) = n / \theta - n \bar{y} ## H(\theta) = U'(\theta) = l''(\theta) = - n / \theta^2 ## I_o(\theta) = - H(\theta) = n / \theta^2 ## I^{-1}_o(\theta) = - H(\theta) = \theta^2/n ## I^{-1/2}_o(\theta) = - H(\theta) = \theta/\sqrt{n} ## Vamos obter de diferentes formas ## i) solução analítica (única necessária nesta caso) ## Demais apenas para ilustração dos métodos: ## ii) solução de equação ## - Método sem uso de derivadas: Brent ## - Método com uso de derivadas: Newton-Raphson ## iii) maximização de função de verossimilhança ## solução analítica (theta.est <- 1/mean(dados)) (lv.theta.est <- lvExp(theta.est, am)) ## estimação numérica ## (somente para ilustração neste caso já que tem-se a analítica) ## (i) via solução de equação (função escore) U(\theta) = l'(\theta) = 0 ## l(\theta) = n log(\theta) - \theta n \bar{y} ## U(\theta) = l'(\theta) = n / \theta - n \bar{y} ## (a) Método de Brent (sem uso de derivadas) UExp <- function(theta, amostra){ return(with(amostra, n*(1/theta - media))) } uniroot(UExp, lower=0, upper=1, amostra=am)$root 1/mean(dados) ## (b) Método de Newton Rapson ## \theta^{r+1} = \theta^{r} - U(\theta) / U'(\theta) ## = \theta^{r} - U(\theta) / H(\theta) ## = \theta^{r} + U(\theta) / I_o(\theta) ## = \theta^{r} - U(\theta) * I^{-1}_o(\theta) IoExp <- function(theta, amostra){ return(amostra$n/theta^2) } ## ou... invIoExp <- function(theta, amostra){ return(theta^2/amostra$n) } maxit <- 100;thetaNR <- 0.1;iter <- 0;d <- 1 while(d > 1e-12 & iter <= maxit){ thetaNR.new <- thetaNR + UExp(thetaNR, am)/IoExp(thetaNR, am) d <- abs(thetaNR - thetaNR.new) thetaNR <- thetaNR.new ; iter <- iter + 1 } thetaNR maxit <- 100;thetaNR <- 0.1;iter <- 0;d <- 1 while(d > 1e-12 & iter <= maxit){ thetaNR.new <- thetaNR + UExp(thetaNR, am) * invIoExp(thetaNR, am) print(c(iter,thetaNR, thetaNR.new,d,maxit)) {curve(lvExp(x, amostra=am), 0.05, 2); abline(v=thetaNR, lty=2); abline(v=thetaNR.new); Sys.sleep(2)} d <- abs(thetaNR - thetaNR.new) thetaNR <- thetaNR.new ; iter <- iter + 1 } thetaNR ## estimativa numérica theta.est ## estimativa analitica curve(UExp(x, am), -0,1);abline(h=0) curve(UExp(x, am), -1,1) UExp(theta.est, am) ## 0 !!! (IoExp.est <- IoExp(theta.est, amostra=am)) (varExp.est <- invIoExp(theta.est, amostra=am)) (sdExp.est <- sqrt(invIoExp(theta.est, amostra=am))) ## verificação de linearidade da score verossimilhança quadrática ## ao redor do máximo ## -I_{-1/2}_o(\hat{\theta}) U(\theta) \approx I_{-1/2}_o(\hat{\theta}) (\theta - \hat{\theta}) (theta.viz <- theta.est + c(-2, 2) * sdExp.est) par(mfrow=c(1,2), mar=c(3,3,0,0), mgp=c(1.8, 0.8, 0)) curve(lvExp(x, amostra=am), 0.05, 2, col=2) abline(v=theta.viz, lty=2) curve(lvExp(x, amostra=am), 0.9*theta.viz[1], 1.1*theta.viz[2], col=2) abline(v=theta.viz, lty=2) par(par.ori) theta.viz <- seq(theta.viz[1],theta.viz[2], l=101) sdUExp <- function(x, amostra) -sdExp.est * UExp(x, amostra=amostra) plot(seq(-2, 2, l=101), sdUExp(theta.viz, amostra=am), type="l"); abline(0,1, lty=3) ## forma alternativa fa <- function(theta, amostra, est){ -sqrt(invIoExp(est, amostra=amostra)) * UExp(theta, amostra=amostra) } fb <- function(theta, amostra, est){ sqrt(IoExp(est, amostra=amostra)) * (theta - est) } plot(fa(theta.viz, amostra=am, est=theta.est) ~ fb(theta.viz, amostra=am, est=theta.est), type="l") abline(0,1, lty=3) ## ## Reparametrização: ## ## 1. Parametrização alternativa da exponencial ## f(y) = (1/\alpha) \exp{- y/alpha} ## l(\alpha) = -n ( \log(\alpha) - \bar{y}/\alpha) ## U(\alpha) = -n (1/\alpha - \bar{y}/\alpha^2) ## H(\alpha) = n (1/\alpha^2 - 2*\bar{y}/\alpha^3) ## = n (\alpha - 2*\bar{y})/\alpha^3 ## I_o(\alpha) = -n (1/\alpha^2 - 2*\bar{y}/\alpha^3) ## = -n (\alpha - 2*\bar{y})/\alpha^3 ## I^{-1}_o(\alpha) = \alpha^3/(n (\alpha - 2*\bar{y})) lvExpa <- function(alpha, amostra) { ## amostra = list(n, media) llik <- with(amostra, -n*(log(alpha) + media/alpha)) return(llik) } curve(lvExpa(x, amostra=am), 0.5, 20) par(mfrow=c(1,2), mar=c(3,3,0,0), mgp=c(1.8, 0.8, 0)) curve(lvExp(x, amostra=am), 0.05, 2, col=2) curve(lvExpa(x, amostra=am), 1/2, 1/0.05, col=2) par(par.ori) ## Obs: ## \alpha = 1/\theta ## limites equivalentes (0.05, 2) --> (0.5, 20) ## invariância da verossimilhança ## \hat{\alpha} = 1/\hat{\theta} UExpa <- function(alpha, amostra){ return(with(amostra, -n*(1/alpha - media/alpha^2))) } ## ou...usando a anterior já programada para parametrização em \theta ## U(\alpha) = U(g(\theta)) * J ## J = d\theta/d\alpha = -1/alpha^2 UExpa1 <- function(alpha, amostra){ theta = 1/alpha return(UExp(theta, amostra = amostra) * (-1/alpha^2)) } ## conferindo que ambas fornecem o mesmo resultado plot(UExpa(2:15, amostra=am), UExpa1(2:15, amostra=am)) abline(0,1) ## I_o(\alpha) = n (\alpha - 2*\bar{y})/\alpha^3 IoExpa <- function(alpha, amostra){ return(with(amostra, -(n/alpha^3) * (alpha-2*media))) } ## ou pela parametrização de theta teriamos que calcular: ## I_o(\alpha) = I_o(g(\theta)) * J^2 + U(g(\theta)) d^2 \theta/d alpha^2 ## ou a inversa invIoExpa <- function(alpha, amostra){ return(with(amostra, -(alpha^3/(n* (alpha-2*media))))) } # examinar diferentes valores inicais no gráfico maxit <- 100;alphaNR <- 0.1;iter <- 0;d <- 1 #maxit <- 100;alphaNR <- 2.5;iter <- 0;d <- 1 #maxit <- 100;alphaNR <- 5;iter <- 0;d <- 1 #maxit <- 100;alphaNR <- 10;iter <- 0;d <- 1 while(d > 1e-12 & iter <= maxit){ alphaNR.new <- alphaNR + UExpa(alphaNR, am)/IoExpa(alphaNR, am) print(c(iter,alphaNR, alphaNR.new,d,maxit)) {curve(lvExpa(x, amostra=am), 0.5, 15); abline(v=alphaNR, lty=2); abline(v=alphaNR.new); Sys.sleep(2)} d <- abs(alphaNR - alphaNR.new) alphaNR <- alphaNR.new ; iter <- iter + 1 } (alpha.est <- mean(dados)) alphaNR 1/thetaNR iter UExpa(alpha.est, am) ## 0 !!! (IoExpa.est <- IoExpa(alpha.est, amostra=am)) (varExpa.est <- invIoExpa(alpha.est, amostra=am)) (sdExpa.est <- sqrt(invIoExpa(alpha.est, amostra=am))) ## verificação de linearidade da score verossimilhança quadrática ## ao redor do máximo ## -I_{-1/2}_o(\hat{\theta}) U(\theta) \approx I_{-1/2}_o(\hat{\theta}) (\theta - \hat{\theta}) (alpha.viz <- alpha.est + c(-2, 2) * sdExpa.est) par(mfrow=c(1,2), mar=c(3,3,0,0), mgp=c(1.8, 0.8, 0)) curve(lvExpa(x, amostra=am), 0.5, 20) abline(v=c(alpha.viz), lty=2) curve(lvExpa(x, amostra=am), 0.9*alpha.viz[1], 1.1*alpha.viz[2]) abline(v=c(alpha.viz), lty=2) par(par.ori) alpha.viz <- seq(alpha.viz[1],alpha.viz[2], l=101) sdUExpa <- function(x, amostra) -sdExpa.est * UExpa(x, amostra=amostra) plot(seq(-2, 2, l=101), sdUExpa(alpha.viz, amostra=am), type="l") abline(0,1, lty=3) ## forma alternativa faa <- function(alpha, amostra, est){ -sqrt(invIoExpa(est, amostra=amostra)) * UExpa(alpha, amostra=amostra) } fba <- function(alpha, amostra, est){ sqrt(IoExpa(est, amostra=amostra)) * (alpha - est) } plot(faa(alpha.viz, amostra=am, est=alpha.est) ~ fba(alpha.viz, amostra=am, est=alpha.est), type="l") abline(0,1, lty=3) ## par(mfrow=c(2,2), mar=c(3,3,0,0), mgp=c(2,1,0)) plot(theta.viz, lvExp(theta.viz, amostra=am), xlab=expression(theta), ylab=expression(L(theta)), type="l") abline(v=range(theta.viz), lty=2) plot(alpha.viz, lvExpa(alpha.viz, amostra=am), xlab=expression(alpha), ylab=expression(L(alpha)), type="l") abline(v=range(alpha.viz), lty=2) plot(seq(-2, 2, l=101), sdUExp(theta.viz, amostra=am), type="l", xlab=expression(z[theta]), ylab=expression(z[U(theta)])) plot(seq(-2, 2, l=101), sdUExpa(alpha.viz, amostra=am), type="l", xlab=expression(z[alpha]), ylab=expression(z[U(alpha)])) par(par.ori) ## ## 2. Parametrizações ## \psi = \log{\theta} ## \phi = \log{\alpha} ## para espaço paramétrico nos reais ## ## l(\psi) = n[log(\theta)-\theta\bar{y}]= n[psi-exp{psi}\bar{y}] ## U(\psi) = l'(g(\theta))*J = n[1-exp{psi}\bar{y}] ## H(\psi) = -n exp{\psi} \bar{y} ## I_o(\psi) = - H(\theta) = n exp{\psi} \bar{y} ## I^{-1}_o(\psi) = exp{-\psi}/(n \bar{y}) lvExp <- function(par, amostra, logpar=FALSE) { ## amostra = list(n, media) if(logpar) theta <- exp(par) else theta <- par llik <- with(amostra, n*(log(theta) - theta * media)) return(llik) } curve(lvExp(x,amostra=am, logpar=T), log(0.05), log(2)) lvExpa <- function(par, amostra, logpar=FALSE) { ## amostra = list(n, media) if(logpar) alpha <- exp(par) else alpha <- par llik <- with(amostra, -n*(log(alpha) + media/alpha)) return(llik) } curve(lvExpa(x, amostra=am, logpar=T), log(1/2), log(1/0.05)) par(mfrow=c(2,2), mgp=c(1.8, 0.8, 0, 0), mar=c(3,3,0,0)) curve(lvExp(x,amostra=am, logpar=F), 0.05, 2, xlab=expression(theta), ylab=expression(l(theta))) curve(lvExpa(x, amostra=am, logpar=F), 1/2, 1/0.05, xlab=expression(alpha), ylab=expression(l(alpha))) curve(lvExp(x,amostra=am, logpar=T), log(0.05), log(2), xlab=expression(psi==log(theta)), ylab=expression(l(psi))) curve(lvExpa(x, amostra=am, logpar=T), log(1/2), log(1/0.05), xlab=expression(phi==log(alpha)), ylab=expression(l(phi))) par(par.ori)