Máxima Verossimilhança
Este é um exercício sobre Máxima Verossimilhança, utilizando os dados de presença e ausência.
Objetivos:
i) Obter estimativa de parâmetros psi e p com base em método de máxima verosimilhança;
ii) Representar graficamente a superfície de verosimilhança.
1) Construa uma base de dados fictícia. Dicas de comandos: sample
, rbinom
,rbind
.
- Esta base de dados pode ser compreendida como uma matriz [i,j], com número i de linhas e j de colunas.
- Como tratam-se de dados de presença e ausência, vamos representar as ausências por zeros, e as presenças por um.
2) Pensar em um modelo simples de ocorrência com dois parâmetros:
- psi é a probabilidade de ocorrência, relativa aos locais (LINHAS);
- p é a detecção (condicionada à ocorrência), relativa ao momento da amostragem (COLUNAS);
- Declarar dois objetos, um para representar o número de linhas da matriz e outro para o número de colunas. Dicas de comandos:
dim, ncol, nrow
.
3) Declarar dois vetores para possíveis valores de psi e p. Dicas de comandos: seq
.
- Lembrem-se de que as probabilidades psi e p variam entre 0 e 1.
4) Declarar uma matriz para guardar os valores de verossimilhança para cada combinação de psi e p. Dicas de comandos: matrix, rep
.
- Esta será uma matriz em que o número de células será equivalente ao número de interações entre os valores criados para psi e para p. Isto é, o número de valores de psi multiplicado ao número de valores de p equivale ao número de células.
5) Calcular valores de verosimilhança para cada par de combinação de psi e p. Dicas de comandos: for, if else
.
- 5.1) um
for
para passar por psi
- 5.2) um
for
(dentro do anterior) para passar por p
- 5.3) crie índices para chamar os valores de psi e p dentro dos loops
for
- 5.4) crie um vetor para guardar os valores de verossimilhança das linhas
- 5.5) mais um
for
para calcular os componentes de verossimilhança correspondentes a cada linha (em cada combinação de psi e p) - <color blue>aqui é que entram as formulinhas</color>
- 5.6) use
if else
. Lembrem-se de que linhas em que tudo é zero levam a uma formulinha diferente das linhas que têm pelo menos uma observação
- 5.7) obtenha o valor de verossimilhança dos dados para a combinação corrente de psi e p com base no resultado das formulinhas para cada linha (não esqueça dos logaritmos)
6) Selecione, na matriz em que você guardou os dados de verossimilhança, os valores dos parâmetros psi e p que maximizam a ocorrência dos dados.Dicas de comandos: which
.
7) Faça uma representação gráfica dos valores de verossimilhança encontrados para cada combinação de psi e p. Dicas de comandos: contour
.
8) Adicione no gráfico o ponto de coordenadas (psi,p), correspondentes aos valores dos parâmetros psi e p que maximizam a ocorrência dos dados. Dicas de comandos: abline
.
Uma solução
################################# ## criando os dados #### ## com parâmetros conhecidos #### ################################# ## Alexandre Adalardo de Oliveira ## 21 de outubro de 2013 # processo de interesse: ocorrência # prob de ocorrenciam psi=0.7 ocorre<-sample(c(1,0),10, prob=c(0.7, 0.3), replace=TRUE) #ou #ocorre<-rbinom(10,1, 0.7) ########################## ## processo de observação ######################### # prob de detecção p=0.4 ######################## # 10 sessões de amostra em cada localidade observa= sample(c(1,0), 100, prob= c(0.4, 0.6), replace=TRUE) detecta<-matrix(observa, ncol=10) ## conectando os dois processos detecta[ocorre==0, ]<-0 amostra<-detecta ############################ ## nome das colunas e linhas rownames(amostra)<-paste("local", 1:10, "") colnames(amostra)<-paste("tempo", 1:10, "") ############################### ## Calculando as probabilidade ## para diferentes parâmetros ############################### nocor<-apply(amostra, 1, sum) ## ocorrência em cada localidade nt<- dim(amostra)[2] ## número de visitas a cada localidade ## sequencia de valores de psi e p seqpsi<-seqp<-seq(0.05,0.955, 0.05) ## criando uma matrix para guardar os resultados resulta<-matrix(NA, ncol=length(seqp), nrow=length(seqpsi), dimnames=list(seqpsi,seqp)) ## um função para o cálculo de verossimilhança para valores de p e psi verosim<-function(psi, p, nocor,nt) { n0 = nocor==0 res<-psi * p ^ nocor * (1-p)^(nt-nocor) res[n0]<-psi * (1-p)^nt + (1-psi) return(res) } ############################################################# ## calculando a Verossimilhança para a sequência de valores ## de p e psi ############################################################# for (i in 1:(length(seqpsi))) { psi<-seqpsi[i] for (j in 1:(length(seqp))) { p<-seqp[j] resulta[i,j]<- -1*sum(log(verosim(psi, p, nocor, nt))) } } ####################################################### ## encontrando a estimativa de Maxima Verossimilhança ## ou minima -LogVerossimilhança ####################################################### indmin<-arrayInd(which.min(resulta), dim(resulta)) psiEst<-dimnames(resulta)[[1]][[indmin[1]]] pEst<-dimnames(resulta)[[2]][[indmin[2]]] psiEst; pEst ###################################################### ## O gráfico da superfície de Verossimilhança contour(resulta, xlab="Ocorrência (psi)", ylab="Detecção (p)", nlevels=20) abline(v=as.numeric(psiEst), col="red", lty=2) abline(h=as.numeric(pEst), col="red", lty=2) ######################################################## # Um conjunto de dados que não conhecemos os parâmetros ######################################################## ## dados de Ploeophana em 10 árvores #################################### ploe<-rbind(c( 0, 0, 0, 0, 0, 0), c( 0, 0, 0, 0, 0, 1), c( 0, 0, 0, 0, 0, 0), c( 0, 0, 0, 0, 0, 0), c( 0, 0, 0, 0, 0, 1), c( 0, 1, 1, 0, 0, 0), c( 0, 1, 1, 1, 0, 0), c( 0, 0, 0, 0, 0, 0), c( 0, 0, 0, 1, 0, 0), c( 1, 0, 1, 0, 1, 0), c( 1, 1, 1, 1, 1, 1)) ##################################################################### # calcule os parâmetros e faça o gráfico da função de verossimilhança #####################################################################