require(gbs)
require(MASS)
require(gamlss.data)
source("BBS-functions.R")

# Nesse código analisamos a funções de verossimilhança de amostras analisadas

# Amostra
# xp = runoff
data(glass)
xp = glass[,1]

# Função penalizada pela priori de Jeffreys modificada
pl = function(x, theta){
  if(theta[1]<0||theta[2]<0){
    f=Inf
    return(f)
  }else{
    a = theta[1]; b = theta[2]; g = theta[3]  
    w = dnorm(g)/pnorm(-g)
    t = (sqrt(x/b) - sqrt(b/x))/a
    f = sum(-log(a*pnorm(-g)*sqrt(32*b*pi))
             -1.5*log(x) + log(x+b) - .5*(abs(t)+g)^2) +
      .5*log( (g-w)*w*(3+g*(g-w))/2 + 1 ) - .5*log(1 + a^2)
    return(f)
  }
}

# Função de verossimilhança original, não penalizada
npl = function(x, theta){
  if(theta[1]<0||theta[2]<0){
    f=Inf
    return(f)
  }else{
    a = theta[1]; b = theta[2]; g = theta[3]  
    w = dnorm(g)/pnorm(-g)
    t = (sqrt(x/b) - sqrt(b/x))/a
    f = sum(-log(a*pnorm(-g)*sqrt(32*b*pi))
             -1.5*log(x) + log(x+b) - .5*(abs(t)+g)^2)
    return(f)
  }
}

# Função de veross. considerando beta fixo no valor verdadeiro
fnpl.ag = function(a,g){
  n = length(a)
  L = rep(0,n)
  for(i in 1:n) L[i] = npl(xp,c(a[i],1.59,g[i]))
  return(L)
} 
# Função de veross. considerando beta fixo no valor verdadeiro
fpl.ag = function(a,g){
  n = length(a)
  L = rep(0,n)
  for(i in 1:n) L[i] = pl(xp,c(a[i],1.59,g[i]))
  return(L)
} 


# combinações dos parâmetros
a = seq(.2,4,length=100)
g = seq(0,20,length=100)

require(lattice)

ag = expand.grid(a,g)

fnpl = fnpl.ag(ag[,1],ag[,2])
wireframe(fnpl~ag[,1]*ag[,2],xlab=expression(alpha),ylab=expression(gamma),
          scales = list(arrows = FALSE),drape=FALSE)
fnpl = outer(a,g,fnpl.ag)
contour(a,g,fnpl,nlevels=200,xlab=expression(alpha),ylab=expression(gamma))
# estimativas obtidas durante o processo de estimação com EMV ordinário
points(c(.27,.33,.5,.67,1.09,1.51,2.79,4.22),
       c(.35,.66,2.18,3.57,6.59,9.38,17.42,26.2),col=1,pch=16)
x0 = c(.27,.33,.5,.67,1.09,1.51,2.79)  # alfa
y0 = c(.35,.66,2.18,3.57,6.59,9.38,17.42)  # gama
x1 = c(.33,.5,.67,1.09,1.51,2.79,4.22)   # alfa
y1 = c(.66,2.18,3.57,6.59,9.38,17.42,26.2)  # gama
# caminho das estimativas
arrows(x0,y0,x1,y1,col=1,lwd=1.5,lty=2)
title(sub="(a)")

fpl = fpl.ag(ag[,1],ag[,2])
wireframe(fpl~ag[,1]*ag[,2],xlab=expression(alpha),ylab=expression(gamma),
          scales = list(arrows = FALSE),drape=FALSE)
fpl = outer(a,g,fpl.ag)
contour(a,g,fpl,nlevels=200,xlab=expression(alpha),ylab=expression(gamma))
title(sub="(b)")

# Gerando o gráfico para o Latex
postscript(file="glass_logliks.eps",horizontal = FALSE, width = 10.0, height = 8.0,
           paper="special")
par(mfrow=c(1,2))
contour(a,g,fnpl,nlevels=200,xlab=expression(alpha),ylab=expression(gamma))
points(c(.27,.33,.5,.67,1.09,1.51,2.79,4.22),
       c(.35,.66,2.18,3.57,6.59,9.38,17.42,26.2),col=1,pch=16)
arrows(x0,y0,x1,y1,col=1,lwd=1.5,lty=2); title(sub="(a)")
fpl = outer(a,g,fpl.ag)
contour(a,g,fpl,nlevels=200,xlab=expression(alpha),ylab=expression(gamma))
title(sub="(b)")
dev.off()


