# # CREDIT TO: 
# https://radfordneal.wordpress.com/2008/08/09/inconsistent-maximum-likelihood-estimation-an-ordinary-example/ 

rm(list=ls())

# GENERATE DATA FROM THE MODEL WITH A GIVEN PARAMETER VALUE.  Arguments are
# the parameter value, t, and the number of data points to generate, m.

gen <- function (t,n)
{
  m <- rep(0,n)
  s <- rep(1,n)
  
  w <- runif(n) < 1/2
  m[w] <- t
  s[w] <- exp(-1/t^2)

  rnorm(n,m,s)
}

Density <- function (x,t)
{
  ll1 <- 0.5*dnorm(x,0,1)

  ll2 <- 0.5*dnorm(x,t,exp(-1/t^2))
#   ll2[x==t] <- -0.5*log(2*pi) + 1/t^2 + log(0.5)

  ll1 + ll2 
}

# COMPUTE LOG DENSITIES OF DATA VALUES FOR A GIVEN PARAMETER VALUE.  Arguments
# are a vector of data values, x, and the parameter value, t.  
#
# Special care is taken to compute the log of extreme density values without 
# overflow.  Data points exactly equal to the parameter value are treated
# specially, since for them the log of the exponential part of the density 
# vanishes, producing results that don't overflow even when the standard
# deviation, exp(-1/t^2), overflows.


log.density <- function (x,t)
{
  ll1 <- dnorm(x,0,1,log=TRUE) + log(0.5)

  ll2 <- dnorm(x,t,exp(-1/t^2),log=TRUE) + log(0.5)
  ll2[x==t] <- -0.5*log(2*pi) + 1/t^2 + log(0.5)

  ll <- rep(NA,length(x))
  for (i in 1:length(x))
  { ll[i] <- add.logs(ll1[i],ll2[i])
  }

  ll
}


# COMPUTE THE LOG LIKELIHOOD GIVEN A DATA VECTOR AND PARAMETER VALUE.  Arguments
# are the vector of data values, x, and the parameter value, t.

log.lik <- function (x,t)
{
  sum(log.density(x,t))
}


# PLOT THE DENSITY FUNCTION FOR A GIVEN PARAMETER VALUE.  Arguments are the
# parameter value, t, and the grid of data values at which to evaluate the
# density.  The grid defaults to 200 points from -3 to 5 (or the parameter
# value, if greater), plus the value of the parameter, where there may be a
# narrow peak in the density that would be missed by the grid.  The vertical
# scale is scaled by a power of ten close to the maximum density, since the
# density values might otherwise overflow.

plot.density <- function (t, grid=sort(c(t,seq(-3,max(5,t),length=200))))
{
  ld <- log.density(grid,t)
  max.ld <- round(max(ld)/log(10))
  ld <- ld - max.ld*log(10)

  plot(grid,exp(ld),type="l",
    xlab="data value",
    ylab=paste("density (x 10^",max.ld,")",sep=""))

  title(paste("Density with parameter",t))
}


# PLOT THE LIKELIHOOD FUNCTION FOR GIVEN DATA, AND FIND MLE.  The argument 
# is a vector of data values, x.  Additional arguments may be given, and are
# passed on the the plot funcition.  The likelihood is plotted for a grid of
# parameter values from 0.01 to 3, in steps of 0.01, plus all positive
# data values less than 3.  An approximation to the MLE is found by taking 
# the maximum over values on this grid.  Data values in (-3,3) are shown at 
# the bottom of the plot. The vertical scale is scaled by a power of ten close 
# to the maximum likelihood, since likelihood values might otherwise overflow.

plot.lik <- function (x,...)
{
  grid <- sort(c(x[x>0 & x<3],seq(0.01,3,by=0.01)))

  ll <- rep(NA,length(grid))
  for (i in 1:length(grid))
  { ll[i] <- log.lik(x,grid[i])
  }

  mlv <- max(ll)
  mle <- grid[ll==mlv][1]
#   max.ll <- round(mlv/log(10))
#   max.ll <- round(log(mlv)/log(10));

#   ll <- ll - max.ll*log(10)
  lik <- exp(ll);
#   lik <- ll;

#   plot(c(-3,3),c(0,max(lik)),type="n",
#   plot(c(-3,3),c(-100,mlv),type="n",
#   plot(c(-3,3),range(ll),type="n",
  plot(c(0,3),range(ll),type="n",
#     xlab="data / parameter",
    xlab="data",
#     ylab=paste("Log-likelihood/10^",max.ll,")",sep=""),
    ylab=paste("Log-likelihood",sep=""),
    ...)
#   points(x,rep(min(ll),length(x)),lty=19)
  lines(grid,ll)

  title(paste("Log-Likelihood -",length(x),"data points,",
              "MLE approximately",round(mle,4)))
}


# ADD NUMBERS REPRESENTED BY LOGS.  Computes log(exp(a)+exp(b)) in a way
# that avoids overflow/underflow problems.

add.logs <- function (a,b)
{ if (a>b)
  { a + log(1+exp(b-a))
  }
  else
  { b + log(1+exp(a-b))
  }
}




fce <- function(x) min(x[x>0])

x <- gen(2,100)

fce(x[1:10])
fce(x[1:30])
fce(x)


# # Plotting observed likelihood
plot.lik(x)

# # Plotting expected likelihood 
theta <- 2; # true value of the parameter
fce <- function(x, t) log.density(x,t) * Density(x,t=theta);

tgrid <- seq(0.1, 4, by=0.01);
N <- length(tgrid);
exp.loglik <- numeric(N)

for(i in 1:N){
	exp.loglik[i] <- integrate(fce, lower=-5, upper=7, t=tgrid[i])$value;
}

plot(tgrid, exp.loglik, type="l")
abline(v=2, lty="dotted");

# The point at which the expected likelihood is maximised 
tgrid[which.max(exp.loglik)]