#========================================================================================== # Aula 18 da disciplina ce223 (20/05/2011) # Função de verossimilhança # Professor Walmes M. Zeviani # www.leg.ufpr.br/ce223 #========================================================================================== #------------------------------------------------------------------------------------------ # caso Poisson ll <- function(lam, x){ # função de verossimilhança da Poisson ## lam é o valor escalar do parâmetro lambda da distribuição Poisson ## x é o vetor de dados da amostra ## o valor retornado é a verossimilhança em lambda sum(x*log(lam)-lam-log(factorial(x))) } ll(2, c(1,2,3)) # escalar e vetor: funciona ll(c(1,2), c(1,2,3)) # vetor e vetor: não funciona #------------------------------------------------------------------------------------------ # uma amostra aleatória e o gráfico da verossilhança amostra <- rpois(20, lambda=7) seq.lam <- seq(0, 15, l=50) seq.ll <- sapply(seq.lam, ll, x=amostra) plot(seq.ll~seq.lam, type="l") #------------------------------------------------------------------------------------------ # encontrar o lambda que maximiza op <- optim(8, ll, x=amostra, control=list(fnscale=-1)) str(op) abline(v=op$par, h=op$value, col="red") #------------------------------------------------------------------------------------------ # função de verossimilhança modificada ll.m <- function(lam, x){ # função de verossimilhança da Poisson ## lam é o valor escalar do parâmetro lambda da distribuição Poisson ## x é o vetor de dados da amostra ## o valor retornado é a verossimilhança em lambda sum(dpois(x, lambda=lam, log=TRUE)) } op.m <- optim(8, ll.m, x=amostra, control=list(fnscale=-1)) str(op.m) abline(v=op$par, h=op$value, col="blue") #------------------------------------------------------------------------------------------ # função de verossimilhança para a beta ll.beta <- function(pars, x){ sh1 <- pars[1]; sh2 <- pars[2] sum(dbeta(x, shape1=sh1, shape2=sh2, log=TRUE)) } amostra.beta <- rbeta(23, 2, 3) amostra.beta op.beta <- optim(c(1.5, 2.5), ll.beta, x=amostra.beta, control=list(fnscale=-1)) op.beta seq.sh1 <- seq(1, 3, l=50) seq.sh2 <- seq(1, 4, l=50) seq.sh <- expand.grid(sh1=seq.sh1, sh2=seq.sh2) seq.ll <- apply(seq.sh, 1, function(x){ ll.beta(x, amostra.beta)}) str(seq.ll) str(seq.sh) require(lattice) wireframe(seq.ll~sh1+sh2, seq.sh, drape=TRUE) levelplot(seq.ll~sh1+sh2, seq.sh) trellis.focus("panel", 1, 1, highlight=FALSE) panel.abline(v=op.beta$par[1], h=op.beta$par[2]) #panel.identify() trellis.unfocus() #------------------------------------------------------------------------------------------