In RobinHankin/hyper2: The Hyperdirichlet Distribution, Mark 2

knitr::opts_chunk$set(echo = TRUE)
library("hyper2")

knitr::include_graphics(system.file("help/figures/hyper2.png", package = "hyper2"))

most of the time is taken calculating the calclikesushi and maxlikechunks. With a sample size of 50 in chunk tryH, and 11 points on the graph, it takes about half an hour to run without cache.

To cite the hyper2 package in publications, please use @hankin2017_rmd.

This short document applies hyper3 formalism to a dataset obtained from the preflib library [@MaWa13a] in which ten different types of sushi were placed in order of preference by 5000 judges. Three of the sushis (maguro, toro, tekka-maki) were made with tuna, and it is reasonable to view these three types as being exchangeable in much the same way as the red and blue buses. To proceed, we need to define equivalence classes of sushi---one for tuna-based sushis and one for everything else. Of course, one can also imagine the opposite phenomenon whereby the judge decides he has had enough tuna for the moment and this disincentivises further tuna choices. This mechanism would be modelled by having $\lambda<1$.

a <- as.matrix(read.table("sushi.txt", skip=22))
set.seed(0)
wanted <- sample(seq_len(nrow(a)),50,replace=TRUE,prob=a[,1])
a <- a[wanted,-1]
colnames(a) <- c("first","second","third","fourth","fifth","sixth","seventh","eighth","ninth","tenth")

sushi_types <- c("ebi", "anago", "maguro", "ika", "uni", "sake", "tamago", "toro", "tekka", "kappa")
a[] <- sushi_types[a]
noquote(a)
sushi_eq_classes <-
    c("ebi"  = 1, "anago" = 1, "maguro" = 2, "ika" = 1, "uni" = 1,
      "sake" = 1, "tamago" = 1, "toro" = 2, "tekka" = 2,
      "kappa" = 1)

make_sushi_H <- function(a, lambda){
    H <- hyper3()

    for(i in seq_len(nrow(a))){
        H <- H + cheering3(a[i,], e=sushi_eq_classes, h=c(1,lambda))
    }
    return(H)
}

H1 <- make_sushi_H(a, lambda=1)
H2 <- make_sushi_H(a, lambda=2)

mH1 <- maxp(H1,give=1)
mH2 <- maxp(H2,give=1)

mH1
mH2

f <- function(lambda){
   H <- make_sushi_H(a, lambda = lambda)
   maxp(H,give=1)$`Log-likelihood`
}

s <- exp(seq(from = log(0.9), to = log(6), len = 11))
L <- sapply(s,f)

(jj <- optimize(f, 2:3, maximum=TRUE))

plot(log(s), L-jj$objective, type='b')

sushi <- make_sushi_H(a, lambda = jj$maximum)

(sushi_maxp <- maxp(sushi))

Package dataset {-}

Following lines create sushi.rda, residing in the data/ directory of the package.

sushi_table <- a
save(sushi_table,  sushi, sushi_maxp, sushi_eq_classes, file="sushi.rda")

References {-}

RobinHankin/hyper2 documentation built on April 13, 2025, 9:33 a.m.