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
#####################################################################