#### Exemplos de como construir verossimilhanças ### ## Exemplo 1 - X ~ N(mu,1) ## 1) Simular do modelo x <- rnorm(100,mean=10,sd=1) ## 2) Construir a verossimilhança ou log.verossimilhança log.vero <- function(mu,x){ log.vero <- sum(dnorm(x,mean=mu,sd=1,log=TRUE)) return(log.vero) } ## 3) Nos casos onde é possível podemos construir um grafico # 3.1) Primeiro construimos um gride para avaliar a função mu.grid <- seq(0,20,by=0.01) # 3.2) Avaliamos a função em todos os pontos do gride valore.log.vero = apply(as.matrix(mu.grid),1,log.vero,x=x) # 3.3) Desenhando o grafico plot(mu.grid,valore.log.vero,type="l",xlab=expression(mu),ylab="Log-Verossimilhança",main="Log-Verossimilhança \n Média da Normal") ## 4) Encontrar o Máximo - No caso geral pode-se fazer usando um algoritmo numerico. estimativas = optim(c(0),log.vero,x=x,control=(list(fnscale=-1)),hessian=TRUE) #### No caso uniparametrico não tem como construir intervalo pela optim. ## Exemplo 1 - X ~ N(mu,sigma^2) ## 1) Simular do modelo x <- rnorm(100,mean=10,sd=10) ## 2) Construir a verossimilhança ou log.verossimilhança log.vero <- function(parametros,x){ mu = parametros[1] sigma = parametros[2] log.vero <- sum(dnorm(x,mean=mu,sd=sigma,log=TRUE)) return(log.vero) } ## 2) Usando um algoritmo numerico estimativas = optim(c(1,1),log.vero,x=x,control=(list(fnscale=-1)),hessian=TRUE) estimativas ## 3) Construindo intervalos de confiança baseado na Io Io <- -estimativas$hessian erro <- qnorm(0.975)*diag(solve(Io)) intervalo.mu = c(estimativas$par[1] - erro[1], estimativas$par[1], estimativas$par[1] + erro[1]) intervalo.mu intervalo.sigma = c(estimativas$par[2] - erro[2], estimativas$par[2], estimativas$par[2] + erro[2]) intervalo.sigma ## 4) Comparando com o resultado analitico media = mean(x) desvio <- sqrt(sum(((x - mean(x))^2)/length(x))) n <- length(x) Io.ana <- matrix(c(n/desvio^2,0,0,2*n/desvio^2), ncol=2,nrow=2) erro.ana = diag(qnorm(0.975)*solve(Io.ana)) erro.ana intervalo.media = c(media - erro.ana[1],media,media+erro.ana[1]) intervalo.desvio = c(desvio - erro.ana[2],desvio, desvio +erro.ana[2]) intervalo.media intervalo.desvio ## 3) X ~ G(a,b) ## 1) Simular do modelo x <- rgamma(1000,shape=1,scale=1/10) ## 2) Construindo a log-verossimilhança log.vero <- function(parametros,x){ a <- parametros[1] b <- parametros[2] log.vero <- sum(dgamma(x,shape=a,scale=1/b,log=TRUE)) return(log.vero) } ## 3) Encontrando o maximo numericamente estimativas = optim(c(1,1),log.vero,x=x,control=(list(fnscale=-1)),hessian=TRUE) estimativas ## 4) Construindo intervalos de confiança baseado na Io Io <- -estimativas$hessian erro <- qnorm(0.975)*diag(solve(Io)) intervalo.a = c(estimativas$par[1] - erro[1], estimativas$par[1], estimativas$par[1] + erro[1]) intervalo.a intervalo.b = c(estimativas$par[2] - erro[2], estimativas$par[2], estimativas$par[2] + erro[2]) intervalo.b