# STA 410, ASSIGNMENT 2, SPRING 2004.


# PLOT DATA AS POINTS ON A CIRCLE.

plot.circle <- function (x)
{
  g <- seq(0,2*pi,length=100)

  plot(cos(g),sin(g),type="l")
  points(cos(x),sin(x),pch=20)
}


# FIND VALUES OF NUMBERS IN A VECTOR MODULO 2PI.

mod2pi <- function (x)
{
  f <- x/(2*pi)
  f <- f-floor(f)
  f*(2*pi)
}


# FIND THE PROBABILITY DENSITY OF POINTS UNDER THE WRAPPED NORMAL DISTRIBUTION.
# The arguments are the the vector of points, the mean of the wrapped normal,
# and the log of the standard deviation of the wrapped normal.  The points and
# the mean are moved to the interval [0,2pi) if necessary.
#
# The density is computed by adding the density of the point under the normal
# distribution with means of mu+2*pi*i for i an integer from -t to t, with t
# chosen so that 2*pi*i will span a range of at least 10 standard deviations.

dwrpnorm <- function (x, mu, rho)
{
  x <- mod2pi(x)
  mu <- mod2pi(mu)
  sd <- exp(rho)

  t <- ceiling(10*sd/(2*pi))
  d <- rep(0,length(x))

  for (i in (-t):t)
  { d <- d + dnorm (x, mu+2*pi*i, sd)
  }

  d
}


# LOG LIKELIHOOD FUNCTION FOR WRAPPED NORMAL MODEL.  Arguments are the vector
# of data points, and the mean and log standard deviation of the distribution.

lglik <- function (x, mu, rho)
{ 
  sum (log (dwrpnorm (x, mu, rho)))
}


# DERIVATIVE OF LOG LIKELIHOOD WITH RESPECT TO MU.

lglik.dmu <- function (x, mu, rho)
{ 
  x <- mod2pi(x)
  mu <- mod2pi(mu)
  sd <- exp(rho)

  t <- ceiling(10*sd/(2*pi))
  s <- 0

  for (i in (-t):t)
  { s <- s + dnorm (x, mu+2*pi*i, sd) * (x-mu-2*pi*i) / sd^2
  }

  sum (s / dwrpnorm (x, mu, rho))
}


# DERIVATIVE OF LOG LIKELIHOOD WITH RESPECT TO RHO.

lglik.drho <- function (x, mu, rho)
{ 
  x <- mod2pi(x)
  mu <- mod2pi(mu)
  sd <- exp(rho)

  t <- ceiling(10*sd/(2*pi))
  s <- 0

  for (i in (-t):t)
  { s <- s + dnorm (x, mu+2*pi*i, sd) * ((x-mu-2*pi*i)^2/sd^2-1)
  }

  sum (s / dwrpnorm (x, mu, rho))
}


# SECOND DERIVATIVE OF LOG LIKELIHOOD WITH RESPECT TO MU.

lglik.dmu2 <- function (x, mu, rho)
{ 
  x <- mod2pi(x)
  mu <- mod2pi(mu)
  sd <- exp(rho)

  t <- ceiling(10*sd/(2*pi))
  s <- 0
  ss <- 0

  for (i in (-t):t)
  { di <- dnorm (x, mu+2*pi*i, sd) 
    ti <- (x-mu-2*pi*i) / sd^2
    s <- s + di*ti
    ss <- ss + di*ti^2
  }

  d <- dwrpnorm (x, mu, rho)

  sum (- s^2 / d^2 + (ss - d/sd^2) / d)
}


# SECOND DERIVATIVE OF LOG LIKELIHOOD WITH RESPECT TO RHO.

lglik.drho2 <- function (x, mu, rho)
{ 
  x <- mod2pi(x)
  mu <- mod2pi(mu)
  sd <- exp(rho)

  t <- ceiling(10*sd/(2*pi))
  s <- 0
  ss <- 0

  for (i in (-t):t)
  { di <- dnorm (x, mu+2*pi*i, sd) 
    ti <- (x-mu-2*pi*i)^2/sd^2
    s <- s + di*(ti-1)
    ss <- ss + di*((ti-1)^2-2*ti)
  }

  d <- dwrpnorm (x, mu, rho)

  sum (- s^2 / d^2 + ss / d)
}


# SECOND DERIVATIVE OF LOG LIKELIHOOD WITH RESPECT TO MU AND RHO.

lglik.dmu.drho <- function (x, mu, rho)
{ 
  x <- mod2pi(x)
  mu <- mod2pi(mu)
  sd <- exp(rho)

  t <- ceiling(10*sd/(2*pi))
  s1 <- 0
  s2 <- 0
  ss <- 0

  for (i in (-t):t)
  { di <- dnorm (x, mu+2*pi*i, sd) 
    ui <- (x-mu-2*pi*i) / sd^2
    vi <- ui * (x-mu-2*pi*i)
    s1 <- s1 + di*ui
    s2 <- s2 + di*(vi-1)
    ss <- ss + di * (vi*ui-3*ui)
  }

  d <- dwrpnorm (x, mu, rho)

  sum (- s1*s2 / d^2 + ss / d)
}


# FIND A MAXIMUM LIKELIHOOD ESTIMATE FOR MU, WITH RHO FIXED.  The arguments
# are the data vector, the initial value for mu, the fixed value for rho,
# the number of iterations of Newton-Raphson iteration to do (default 20),
# and a debug flag (default FALSE), which when TRUE causes the parameter
# value and log likelihood to be printed for every iteration.  

# The value returned is a possible maximum likelihood estimate for mu.  This
# value will be a zero of the derivative of the log likelihood, if the
# Newton-Raphson iteration has converged.  No attempt is made here to check
# whether it is actually a maximum (rather than a minimum or inflection point).

mle.mu <- function (x, mu, rho, n=20, debug=FALSE)
{ 
  if (debug)
  { cat(0,mu,lglik(x,mu,rho),"\n")
  }

  for (k in 1:n)
  { 
    # Newton-Raphson iteration, with mu moved back to [0,2p) afterwards.

    mu <- mod2pi (mu - lglik.dmu(x,mu,rho) / lglik.dmu2(x,mu,rho))

    if (debug)
    { cat(k,mu,lglik(x,mu,rho),"\n")
    }
  }

  mu
}


# FIND A JOINT MAXIMUM LIKELIHOOD ESTIMATE FOR MU AND RHO.  The arguments
# are the data vector, the initial values for mu and rho, the number of 
# iterations of Newton-Raphson iteration to do (default 20), and a debug 
# flag (default FALSE), which when TRUE causes the parameter values and log 
# likelihood to be printed for every iteration.  

# The value returned is a list with elements "mu" and "rho" containing the 
# possible maximum likelihood estimates and an element "cov" containing the
# estimated covariance matrix for these maximum likelhood estimates.  The
# square roots of the diagonal elements of the "cov" matrix can be used as 
# standard errors for the estimates.  The values for mu and rho returned will 
# be a point where the derivatives of the log likelihood are zero, if the
# Newton-Raphson iteration has converged.  No attempt is made here to check
# that this is actually a maximum (rather than a minimum or saddle point).

mle.joint <- function (x, mu, rho, n=20, debug=FALSE)
{
  if (debug)
  { cat(0,mu,rho,lglik(x,mu,rho),"\n")
  }

  H <- matrix(0,2,2)

  for (k in 1:n)
  { 
    # Put current parameters into a vector for use in N-R iteration.

    v <- c(mu,rho)

    # Put second derivatives into a matrix (the Hessian).

    H[1,1] <- lglik.dmu2(x,mu,rho)
    H[2,2] <- lglik.drho2(x,mu,rho)
    H[1,2] <- H[2,1] <- lglik.dmu.drho(x,mu,rho)

    # Put first derivatives into a vector (the gradient).

    g <- c (lglik.dmu(x,mu,rho), lglik.drho(x,mu,rho))

    # Do the Newton-Raphson iteration.

    v <- v - solve(H,g)

    # Take parameters out of the resulting vector, and move mu into [0,2pi).

    mu <- mod2pi(v[1])
    rho <- v[2]

    if (debug)
    { cat(k,mu,rho,lglik(x,mu,rho),"\n")
    }
  }

  list (mu=mu, rho=rho, cov=solve(-H))
}