R/LInormVar.R

In LBI: Likelihood Based Inference

LInormVar = function(x, k, conf.level=0.95)
{
  x = x[!is.na(x)]
  n0 = length(x) ; sxx = sum(x^2) ; sx = sum(x)
  if (!is.numeric(x) | sum(is.infinite(x) > 0) | sum(is.nan(x)) > 0 | n0 < 3 | length(unique(x)) == 1) stop("Check the input!")
  m0 = sx/n0
  v0 = sxx/n0 - m0^2
  s0 = sqrt(v0)
  maxLL = -n0*(log(2*pi*v0) + 1)/2

  if (!missing(k)) {
    logk = log(k)
  } else {
    logk = n0/2*log(1 + qf(conf.level, 1, n0 - 2)/(n0 - 2)) # two parameters with one nuisance (mean)
    logk = min(logk, log(2/(1 - conf.level)))
  }

  O2 = function(th) maxLL + (n0*log(2*pi*th) + (sxx - sx^2/n0)/th)/2 - logk
  varLL = uniroot(O2, c(1e-8, v0))$root
  varUL = uniroot(O2, c(v0, 100*v0))$root
  sdLL = sqrt(varLL)
  sdUL = sqrt(varUL)
  Res = cbind(PE = c(s0, v0), LL=c(sdLL, varLL), UL=c(sdUL, varUL))
  rownames(Res) = c("sd", "var")
  attr(Res, "n") = n0
  attr(Res, "k") = exp(logk)
  attr(Res, "log(k)") = logk
  attr(Res, "maxLogLik") = maxLL
  return(Res)
}