#==========================================================================================
# 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()

#------------------------------------------------------------------------------------------