R/foldedt.mle.R

In foldedt: The Folded t Family of Distributions

foldedt.mle <- function(x) {

  fun <- function(para, n, x) {
    mu <- para[1]  ;  s2 <- exp(para[2])  ;  v <- exp(para[3])
    vs2 <- v * s2
    a1 <- ( 1 + (x - mu)^2/vs2 )^(- 0.5 * (v + 1) )
    a2 <- ( 1 + (x + mu)^2/vs2 )^(- 0.5 * (v + 1) )
    - n * lgamma(0.5 * v + 0.5) + n * 0.5 * log(pi * vs2) + n * lgamma(0.5 * v) - sum( log(a1 + a2) )
  }
  n <- length(x)
  mod <- optim( c( mean(x), var(x), log(5) ), fun, n = n, x = x, control = list(maxit = 5000) )
  mod <- optim( mod$par, fun, n = n, x = x, control = list(maxit = 5000) )
  param <- c( mod$par[1], exp(mod$par[2:3]) )
  names(param) <- c("location", "scatter", "nu")
  list( param = param, loglik = -mod$value )
}