shifted_log_logistic.R

x <- c(...) # samples

# for fitdist
suppressWarnings(library(fitdistrplus))
# for llog3 and lnorm3
suppressWarnings(library(FAdist))

print_ks_test <- function (x, dist, args) {
  print(do.call(function (...) {
    ks.test(x, paste0("p", dist), ...)
  }, args))
}

print_list_arg <- function (list) {
  m <- matrix(list, dimnames=list(names(list)))
  colnames(m) <- deparse(substitute(list))
  print(m)
}

estimate_thres <- function (x) {
  lower = 0
  upper = min(x)-1
  start = upper * 0.9
  optim(start, function (shift) {
    abs(mean(log(x - shift)) - median(log(x - shift)))
  }, lower=lower, upper=upper, method="Brent")$par
}

fit_dist <- function (x, dist, ...) {
  cat(paste0("\nFit ", dist, "\n"))

  start = list(...)

  print_list_arg(start)

  print_ks_test(x, dist, start)

  fit <- fitdist(x, dist, "mge", start = start, gof="AD")

  print(fit)

  print_ks_test(x, dist, as.list(fit$estimate))

  svg(paste0(dist, ".svg"))
  suppressWarnings(plot(fit))
  dev.off()
}

thres <- estimate_thres(x)
sdlog <- sd(log(x - thres))
meanlog <- mean(log(x - thres))

fit_dist(x, "lnorm3",
  shape = sdlog,
  scale = meanlog,
  thres = thres
)

fit_dist(x, "llog3",
  shape = sdlog * 0.5513289, # sdlog * sqrt(3) / pi
  scale = meanlog,
  thres = thres
)