######### Figuras Capitulo 2 apostila John Tawn ################################ ## Figura 1 - Gráfico típico de verossimilhança vero <- function(theta){ n*log(theta) - theta*sum(x)} set.seed(100) x <- rexp(100) n <- length(x) theta <- seq(0.5,1.5,by=0.001) l.theta <- vero(theta) plot(theta,l.theta,type="l",xlab=expression(theta), ylab="Log-Verossimilhanca") arrows(x1=1/mean(x),y1=min(l.theta),x0=1/mean(x),y0=max(l.theta)) text(x=1/mean(x),y=min(l.theta),round(1/mean(x),digits=2),pos=1) abline(v=1/mean(x)) abline(h=-100) ## Figura 2 - Gráfico atípicos da verossimilhança ## Distribuição Uniforme x <- c(7.1,2.6,0.5,1.3,9.2,3,0.7,1.1,5,1.3) n <- 10 vero <- function(theta){ maximo <- max(x) ifelse(theta >= maximo, theta^-n,0)} theta <- seq(from=8,to=12,by=0.01) l.theta <- vero(theta) plot(theta,l.theta,type="l",ylab="Verossimilhanca",xlab="Valores de theta") ## Verossimilhança bimodal n <- 12 x <- 1:n y <- c(1,3,5,10,8,5,6,7,5,3,2,1) plot(x, y, ylab = "Verossimilhanca", xlab="Theta") lines(spline(x, y)) abline(h=6) ## Exemplo 2.1 x <- rnorm(10,mean=10) theta <- seq(from=8,to=12,by=0.001) vero.normal <- c() for(i in 1:length(theta)){ vero.normal[i] <- -(n/2)*log(2*pi) - 0.5*sum(x-theta[i])^2} plot(theta,vero.normal,type="l",xlab="Valores de theta",ylab="Log-Verossimilhanca") ## Estimativa de maxima verossimilhanca mean(x) abline(v=mean(x)) theta.l <- mean(x) - 1.96/sqrt(length(x)) theta.u <- mean(x) + 1.96/sqrt(length(x)) abline(v=c(theta.l,theta.u)) ## Grafico da deviance theta <- seq(from=9,to=11,by=0.001) deviance <- c() for(i in 1:length(theta)){ deviance[i] <- sum(x - theta[i])^2 - sum(x - mean(x))^2} plot(theta,deviance,type="l",xlab="Valores de theta", ylab="Deviance") ######################################################################## ## Exemplo 2.3 ## Impacto do tamanho da amostra na verossimilhança n1 <- 10 n2 <- 25 n3 <- 100 d.theta <- function(theta,n){ 2*n *(theta + 2*(log(2) - log(theta) -1))} theta <- seq(from=1,to=3,by=0.01) dev.theta1 <- d.theta(theta,n1) dev.theta2 <- d.theta(theta,n2) dev.theta3 <- d.theta(theta,n3) plot(theta,dev.theta1,type="l") lines(theta,dev.theta2,type="l",col="red") lines(theta,dev.theta3,type="l",col="blue") text(x=1.5,y=2,"n = 10") text(x=2.4,y=2,"n = 25") text(x=1.8,y=2,"n = 100") ##### Aproximação assintótica da log-verossimilhança caso exponencial set.seed(100) x <- rexp(1000) n <- length(x) theta <- seq(0.5,2,by=0.001) log.vero <- function(theta){ n*log(theta) - theta*n*mean(x)} log.vero.apr <- function(theta){ theta.est <- 1/mean(x) l.the <- n*log(theta.est) - theta.est*sum(x) l1.the <- n/theta.est - n*mean(x) l2.the <- -(n/theta.est^2) l.the + (theta - theta.est)*l1.the + 0.5*(theta - theta.est)^2 * l2.the } plot(theta,log.vero(theta),type="l",main="Aproximacao da log-verossimilhanca") lines(theta,log.vero.apr(theta),type="l",col="red") ## Exemplificando a invariância ## Caso exponencial x <- c(2.063,1.883,2.144,2.491,1.677,3.365,0.979,2.549,0.318,1.704,0.926,2.671) n <- length(x) ## Construindo as duas log-verossimilhanca ## Funcao de theta vero1 <- function(theta){ n*log(2) + n*log(theta) + sum(x) - theta*sum(x^2)} ## Funcao de mu vero2 <- function(mu){ n*log(2) - n*log(mu) - sum(x^2)/mu} ## Grade de valores para os parametros theta <- seq(0,1.8,by=0.001) mu <- seq(1.5,15,by=0.001) ## Plotando a verossimilhança par(mfrow=c(1,2)) plot(theta,vero1(theta),type="l",ylab="Verossimilhanca", xlab=expression(theta)) ### Encontrando estimativas pontual e intervalar ## Para theta theta.est <- length(x)/sum(x^2) abline(v=theta.est) V.theta <- theta.est^2/n V.theta IC.theta <- c(theta.est - 1.96*sqrt(V.theta), theta.est + 1.96*sqrt(V.theta)) abline(v=c(IC.theta)) ### Para mu plot(mu,vero2(mu),type="l",ylab="Verossimilhanca", xlab=expression(mu)) mu.est <- sum(x^2)/n abline(v=mu.est) V.mu <- mu.est^2/n IC.mu <- c(mu.est - 1.96*sqrt(V.mu), mu.est + 1.96*sqrt(V.mu)) abline(v=c(IC.mu)) ## Usando o método Delta para encontrar a variancia de mu dado a de theta # Supondo que eu tenho o intervalo para theta mais meu interesse é em mu = theta^-1 # Temos que mu = g(theta) = theta^-1 logo g'(theta) = theta^-2 # Usando resultados assintóticos temos que gl1.theta <- -1/theta.est^2 V.mu.theta <- (gl1.theta^2)*V.theta V.mu.theta IC.mu.theta <- c(1/theta.est - 1.96*sqrt(V.mu.theta), 1/theta.est + 1.96*sqrt(V.mu.theta)) ## Comparando o obtido direto da verossimilhança com o indireto vindo do metodo Delta IC.mu.theta IC.mu ################################################################################ ## Obtendo a de theta dado mu theta = mu^-1 V.theta.mu <- (1/mu.est^2)^2 * V.mu V.theta.mu V.theta ################################################################################ ## Exemplo de teste de razão de verossimilhanças ## Exemplo X ~ N(theta,1) n =100 testar se theta=6 com 5% de significância x <- rnorm(10,mean=6) ## Construindo a Deviance deviance <- function(theta){ length(x)*(mean(x) - theta)^2} teste <- seq(5.5,7,by=0.1) p.valor <- c() for(i in 1:length(teste)){ p.valor[i] <- 1-pchisq(deviance(teste[i]),df=1)} plot(teste,p.valor,type="l") abline(h=.95)