hdi_TJ.R

compute_hpdi_TJ <- function(xs, prob = .9) {
  x_sorted <- sort(xs)
  n <- length(xs)
  
  num_to_keep <- ceiling(prob * n)
  num_to_drop <- n - num_to_keep

  possible_starts <- seq(1, num_to_drop + 1, by = 1)
  # Just count down from the other end
  possible_ends <- rev(seq(from = n, length = num_to_drop + 1, by = -1))

  # Find smallest interval
  span <- x_sorted[possible_ends] - x_sorted[possible_starts]
  edge <- which.min(span)
  edges <- c(possible_starts[edge], possible_ends[edge])

  # My requirement: length of span interval must be same as number to keep.
  # Other methods produce intervals that are 1 longer.
  stopifnot(length(edges[1]:edges[2]) == num_to_keep)
  
  x_sorted[edges]
}

is.inside <- function(x, i){
  x>=i[1] & x<=i[2]
} 

try_seed_fun <- function(s, rf, df){
  
  set.seed(s)
  x <- rf(1e3)
  d <- df(x)
  
  i_hdi <- HDInterval::hdi(x, credMass=0.8)
  
  i_TJ <- compute_hpdi_TJ(x, prob = 0.8)
  
  tibble::tibble(
    # is TJ's interval smaller
    TJ_smaller= unname(diff(i_hdi))>=diff(i_TJ),
    # are all densities inside higher than densities outside?
    HDI_all_densities=min( d[is.inside(x,i_hdi)] ) > max( d[!is.inside(x, i_hdi)] ),
    TJ_all_densities=min( d[is.inside(x,i_TJ)] ) > max( d[!is.inside(x, i_TJ)] )
  )
}

library(tidyverse)

purrr::map_dfr(1:1000, try_seed_fun, rlnorm, dlnorm) %>% 
  count(TJ_smaller, HDI_all_densities, TJ_all_densities)

# # A tibble: 4 x 4
# TJ_smaller HDI_all_densities TJ_all_densities     n
# <lgl>      <lgl>             <lgl>            <int>
# 1 TRUE       FALSE             FALSE              602
# 2 TRUE       FALSE             TRUE               134
# 3 TRUE       TRUE              FALSE              132
# 4 TRUE       TRUE              TRUE               132