# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo:
#
# Aproximacao de Laplace
#
# Animais!
#==============================================================================
L  <- function(theta,k){
	# k=0 => L1
	# k=1 => L2 
	# k=2 => L3
	y   <- c(125,38,34)
	out <- ( y[1]*log(2+theta)
              +y[2]*log(1-theta)
              +(y[3]+k)*log(theta) )
	return(out)
}

#==============================================================================
# Segunda derivada
#==============================================================================
D2  <- function(theta,k){
	y   <- c(125,38,34)
	out <- ( -y[1]/((2+theta)^2)
               -y[2]/((1-theta)^2)
               -(y[3]+k)/(theta^2))
	return(out)
}

#==============================================================================
# Encontrar modas utilizando grade de valores
#==============================================================================
theta.aux <- seq(0.00001,0.99999,0.00001)
m         <- length(theta.aux)

f1        <- rep(NA,m)
for (i in 1:m){
	f1[i] <- L(theta.aux[i],0)
}
theta1 <- theta.aux[f1==max(f1)]
print(theta1)

f2        <- rep(NA,m)
for (i in 1:m){
	f2[i] <- L(theta.aux[i],1)
}
theta2 <- theta.aux[f2==max(f2)]
print(theta2)

f3        <- rep(NA,m)
for (i in 1:m){
	f3[i] <- L(theta.aux[i],2)
}
theta3 <- theta.aux[f3==max(f3)]
print(theta3)

[1] 0,62682
[1] 0,63099
[1] 0,63506

#==============================================================================
# Psi^2
#==============================================================================
var1     <- (-1/D2(theta1,0))
var2     <- (-1/D2(theta2,1))
var3     <- (-1/D2(theta3,2))

#==============================================================================
# Aproximacao da media a posteriori de theta
#==============================================================================
E.theta  <- exp(L(theta2,1)-L(theta1,0))*sqrt(var2/var1)
print(E.theta)

[1] 0,6227461

#==============================================================================
# Aproximacao do segundo momento ordinario a posteriori de theta
#==============================================================================
E.theta2 <- exp(L(theta3,2)-L(theta1,0))*sqrt(var3/var1)
print(E.theta2)

[1] 0,3904036

#==============================================================================
# Variancia e desvio padrao a posteriori (aproximacao)
#==============================================================================
var.theta <- E.theta2 - E.theta^2
print(var.theta)
print(sqrt(var.theta))

[1] 0,00259099
[1] 0,05090177

#==============================================================================
# Utilizando integracao Monte Carlo
#==============================================================================
set.seed(12345)      # Semente
N       <- 10000     # Tamanho da amostra (Monte Carlo)
rtheta  <- runif(N)
py      <- rep(NA,N)
num1    <- rep(NA,N) 
num2    <- rep(NA,N)

for(i in 1:N){
	py[i]   <- exp(L(rtheta[i],0))
}

for(i in 1:N){ 
	num1[i] <- exp(L(rtheta[i],1))
}

for(i in 1:N){ 
	num2[i] <- exp(L(rtheta[i],2))
}

#==============================================================================
# Aproximacao da media a posteriori de theta
#==============================================================================
E.rtheta  <- sum(num1)/sum(py)
print(E.rtheta)

[1] 0,6220711

#==============================================================================
# Aproximacao do segundo momento ordinario a posteriori de theta
#==============================================================================
E.rtheta2 <- sum(num2)/sum(py)
print(E.rtheta2)

[1] 0,3896079

#==============================================================================
# Variancia e desvio padrao a posteriori (aproximacao)
#==============================================================================
var.rtheta <- E.rtheta2 - E.rtheta^2
print(var.rtheta)
print(sqrt(var.rtheta))

[1] 0,002635459
[1] 0,05133673

# rm(list=ls()) 
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================