###
###  MODERN STATISTICAL METHODS (NMST434)
###  4. exercise class
###
# setwd("H:/nmst434")
 
### In K10, one has to add a net directory L:/Statistika/R/library to the path
### where the R packages are looked for:
.libPaths("L:/Statistika/R/library")      ### This is not needed on most private PC's
                                          ### where R and all extension packages are
                                          ### installed on one place.

# # MLE incosistent estimator 

# Data is i.i.d. real, with distribution (1/2)N(0,1) + (1/2)N(theta,exp(-2/theta^2)),
# where theta is a positive real parameter.

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

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


# COMPUTE LOG DENSITIES OF DATA VALUES FOR A GIVEN PARAMETER VALUE.  Arguments
# are a vector of data values, x, and the parameter value, t.  


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 
}


# 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 <- log(dnorm(x,0,1)) + log(0.5)

  ll2 <- log(dnorm(x,t,exp(-1/t^2))) + 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)))
	
	lines(rep(2,2), range(ll), col="blue", lty="dotted")
}


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


# PRODUCE DEMO PLOTS.  Leaves the generated data in the global variable x for
# possible later examination.

demo <- function()
{ 
  par(mfrow=c(2,3))

  plot.density(2.5)
  plot.density(0.6)
  plot.density(0.2)

#   set.seed(1)
#   x <- gen(0.6,100)
  x <<- gen(2,100)
  plot.lik(x[1:10]);
  plot.lik(x[1:30])
  plot.lik(x)
	
	par(mfrow=c(1,1))
}

demo()

fce <- function(x) min(x[x>0])
fce(x[1:10])
fce(x[1:30])
fce(x)

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