####################
### exercise 2.1 ###
####################
rm(list = ls())

# data form iid Cauchy(theta, 1)
y <- c(1.77, -0.23, 2.76, 3.80, 3.47, 56.75, -1.34, 4.24, -2.44, 3.29, 3.71,
       -2.40, 4.53, -0.07, -1.05, -13.87, -2.53, -1.75, 0.27, 43.21)

### 1.a ###
llik <- function(theta, y)
{
  n <- length(y)
  -n*log(pi) - sum(log(1 + (y - theta)^2))
}

theta.grid <- seq(-60, 60, by = 0.001)
llik.out <- sapply(theta.grid, llik, y = y)
plot(theta.grid, llik.out, type = 'l', ylab = "log-likelihood values",
     xlab = expression(theta))
abline(v = theta.grid[which.max(llik.out)], lty = 2)
theta.grid[which.max(llik.out)]

# the strange shape is not due to the grid
llik <- function(theta)
{
  data <- c(1.77, -0.23, 2.76, 3.80, 3.47, 56.75, -1.34, 4.24, -2.44, 3.29, 3.71,
            -2.40, 4.53, -0.07, -1.05, -13.87, -2.53, -1.75, 0.27, 43.21)
  n <- length(data)
  -n*log(pi) - sum(log(1 + (data - theta)^2))
}
LLIK <- Vectorize(llik)
curve(LLIK, -60, 60)
abline(v = theta.grid[which.max(llik.out)], lty = 2)

# score function
llik.theta <- function(theta, y)
{
  2 * sum((y - theta)/(1 + (y - theta)^2))
}

# observed Fisher information
llik.thetatheta <- function(theta, y)
{
  - 2 * sum( (1 - (y - theta)^2) / ((1 + (y - theta)^2)^2) )
}

# optimizer
NewtonRaphson <- function(x, y, epsilon = 1e-6)
{
  err <- 1e3
  while (abs(err) > epsilon)
  {
    h <- llik.theta(x[length(x)], y) / llik.thetatheta(x[length(x)], y)
    x <- c(x, x[length(x)] - h)
    if(is.infinite(x[length(x)])) return(x)
    err <- x[length(x)] - x[length(x) - 1]
  }
  return(x)
}

NewtonRaphson(-11, y)
NewtonRaphson(-1, y)
NewtonRaphson(0, y)
NewtonRaphson(1.5, y)
NewtonRaphson(4, y)
NewtonRaphson(4.7, y)
NewtonRaphson(7, y)
NewtonRaphson(8, y)
NewtonRaphson(38, y)

x.minus11 <- NewtonRaphson(-11, y)
theta.grid <- seq(-60, 60, by = 0.01)
llik.theta.out <- sapply(theta.grid, llik.theta, y = y)
plot(theta.grid, llik.theta.out, type = 'l', ylab = "score function",
     xlab = expression(theta))
abline(h = 0, lty = 2)
for (i in 1:5)
{
  points(x.minus11[i], llik.theta(x.minus11[i], y))  
  invisible(readline())
}

plotSteps <- function(res, nstep = 5)
{
  # make the plot
  theta.grid <- seq(-60, 60, by = 0.01)
  llik.theta.out <- sapply(theta.grid, llik.theta, y = y)
  plot(theta.grid, llik.theta.out, type = 'l', ylab = "score function",
       xlab = expression(theta))
  abline(h = 0, lty = 2)
  # add the points
  for (i in 1:nstep)
  {
    points(res[i], llik.theta(res[i], y))  
    invisible(readline())
  }
}

x.4.7 <- NewtonRaphson(4.7, y)
plotSteps(x.4.7)

x.meany <- NewtonRaphson(mean(y), y)
plotSteps(x.meany, 10)



### 1.b ###
bisection <- function(x.low, x.upp, y, epsilon = 1e-6)
{
  err <- 1e3
  if(llik.theta(x.low, y) * llik.theta(x.upp, y) > 0)
    stop('Initial conditions not satisfied' )
  x <- x.low + (x.upp - x.low)/2
  while (abs(err) > epsilon)
  {
    ifelse(llik.theta(x.low, y)*llik.theta(x[length(x)], y) < 0,
           x.upp <- x[length(x)], x.low <- x[length(x)])
    x <- c(x, x.low + (x.upp - x.low)/2)
    err <- x[length(x)] - x[length(x) - 1]
  }
  return(x)
}

bisection(-1, 1, y)

theta.grid <- seq(-2, 5, by = 0.01)
llik.theta.out <- sapply(theta.grid, llik.theta, y = y)
plot(theta.grid, llik.theta.out, type = 'l', ylab = "score function",
     xlab = expression(theta))
abline(h = 0, lty = 2)

# bisection method fails when g'(a)*g'(b) > 0. E.g.:
bisection(-1, 2, y)

# bisection method fails to find the global minimum when a local minimum is
# closer then the global minimum to the middle point. E.g.:
bisection(-1, 5, y)



### 1.c ###
fixed.point <- function(x, y, alpha = 1, epsilon = 1e-6, max.iter = 500)
{
  iter <- 1
  err <- 1e3
  while (abs(err) > epsilon && iter < max.iter)
  {
    x <- c(x, alpha * llik.theta(x[length(x)], y) + x[length(x)])
    err <- x[length(x)] - x[length(x) - 1]
    iter <- iter + 1
  }
  return(x)
}

fixed.point(-1, y, 1)
fixed.point(-1, y, 0.64)
fixed.point(-1, y, 0.25)
fixed.point(-1, y, 0.64, epsilon = 1e-3)

fixed.point(0, y, 1)
fixed.point(0, y, 0.64)
fixed.point(0, y, 0.25)
fixed.point(3, y, 1)
fixed.point(3, y, 0.25)

fixed.point(-1, y, 0.5)
fixed.point(-1, y, 0.1)
fixed.point(-1, y, 0.01)



### 1.d ###
secant <- function(x, y, epsilon = 1e-6)
{
  err <- 1e3
  if (length(x) < 2) stop("the secant method requires 2 starting points")
  while (abs(err) > epsilon)
  {
    x <- c(x, x[length(x)] - llik.theta(x[length(x)], y) * 
             (x[length(x)] - x[length(x) - 1])/
             (llik.theta(x[length(x)], y) - llik.theta(x[length(x) - 1], y)))
    err <- x[length(x)] - x[length(x) - 1]
  }
  return(x)
}

secant(c(-2, -1), y)
secant(c(-3, 3), y)
secant(c(-11, -13), y)



### 1.e ### 
# time
time.nr <- system.time(x.nr <- NewtonRaphson(4.7, y))
time.bs <- system.time(x.bs <- bisection(-1, 1, y))
time.fp <- system.time(x.fp <- fixed.point(-1, y, 0.25))
time.se <- system.time(x.se <- secant(c(-2, -1), y))
rbind(time.nr, time.bs, time.fp, time.se)

# iterations
rbind(length(x.nr), length(x.bs), length(x.fp), length(x.se))





####################
### exercise 2.2 ###
####################
rm(list = ls())

y <- c(3.91, 4.85, 2.28, 4.06, 3.70, 4.04, 5.46, 3.53, 2.28, 1.96,
       2.53, 3.88, 2.22, 3.47, 4.82, 2.46, 2.99, 2.54, 0.52, 2.50)

### 2.a ###
llik <- function(theta, y)
{
  n <- length(y)
  -n*log(2 * pi) + sum(log(1 - cos(y - theta)))
}

theta.grid <- seq(-pi, pi, length = 10000)
llik.out <- sapply(theta.grid, llik, y = y)
plot(theta.grid, llik.out, type = 'l', ylab = "log-likelihood values",
     xlab = expression(theta))
abline(v = theta.grid[which.max(llik.out)], lty = 2)
theta.grid[which.max(llik.out)]


### 2.b ###
# E_X[X] = sin(theta) + pi
theta.hat.moments <- asin(mean(y)-pi)


### 2.c ###
llik.theta <- function(theta, y)
{
  - sum( sin(y - theta) / (1 - cos(y - theta) ) )
}

# observed Fisher information
llik.thetatheta <- function(theta, y)
{
  - sum( 1 / (1 - cos(y - theta) ) )
}

NewtonRaphson(theta.hat.moments, y)
NewtonRaphson(-2.7, y)
NewtonRaphson(2.7, y)


### 2.d ###
starting.points <- seq(-pi, pi, le = 200)
onlyOpt <- function(x, y)
{
  tmp <-  NewtonRaphson(x, y)
  tmp[length(tmp)]
}
res <- round(sapply(starting.points, onlyOpt, y = y), 2)
val.found <- unique(res)
groups <- c(min(starting.points[res == val.found[1]]),
            max(starting.points[res == val.found[1]]))
for (i in 2:length(val.found))
  groups <- rbind(groups,
                  c(min(starting.points[res == val.found[i]]),
                    max(starting.points[res == val.found[i]]))
                 )
groups


### 2.e ###
theta.grid[which.min(llik.out)]
NewtonRaphson(theta.grid[which.min(llik.out)] - 0.0001, y)
NewtonRaphson(theta.grid[which.min(llik.out)] + 0.0001, y)


# plot score
theta.grid <- seq(-pi, pi, by = 0.01)
llik.theta.out <- sapply(theta.grid, llik.theta, y = y)
plot(theta.grid, llik.theta.out, type = 'l', ylab = "score function",
     xlab = expression(theta))
abline(h = 0, lty = 2)




### for extra exercise 3f ###
f <- function(x) -x^2 + x + 1/4
x.grid <- seq(-1, 2, length = 1000)
out <- sapply(x.grid, f)
plot(x.grid, out, type = 'l')
abline(v = -1/2, lty = 2)
abline(v = 1/2, lty = 2)
abline(v = 3/2, lty = 2)