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

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