In RobinHankin/hyper2: The Hyperdirichlet Distribution, Mark 2

knitr::opts_chunk$set(echo = TRUE)
library(hyper2)
library(stringr)
set.seed(0)

Basic idea is to have $n$ competitors with BT strengths $\beta, \beta x,\ldots \beta x^{n-1}$ with $x\in[0,1]$. The unit sum constraint gives $\beta=\frac{1-x}{1-x^n}$ but this usually cancels out. One of the competitors is the focal competitor with BT strength $\beta x^{r-1}$ for some $r$ with $1\leqslant r < n$. Thus $r=1$ means that the focal competitor is the strongest among the $n$, $r=2$ means he is the second, and so on.

BT <- function(n,x){
    out <- x^(0:(n-1))
    if(x != 1){
        out <- out*(1-x)/(1-x^n)
    } else {
        out <- out/n
    }
    names(out) <- paste0("p",str_pad(1:n,ceiling(log10(n))))
    out
}
BT(4,0.5)
sum(BT(4,0.5))

We will consider $n=9$, $x=0.5$, $r=4$ and make repeated observations:

(jj <- BT(9,0.5))
o <- tabulate(replicate(1e4,which(rrace(jj) == "p4")))
o
plot(o,type='h')

getprob <- function(n, r, x, N=1e4){
    jj <- BT(n, x)
    wanted <- paste0("p", str_pad(r ,ceiling(log10(n))))
    c(r=r, x=x, setNames(tabulate(replicate(N, which(rrace(jj) == wanted)),nbins=n), paste0("p", 1:n)))
}
rbind(
getprob(9,1,0.5,N=100),
getprob(9,2,0.5,N=100),
getprob(9,3,0.5,N=100),
getprob(9,4,0.5,N=100),
getprob(9,5,0.5,N=100),
getprob(9,6,0.5,N=100),
getprob(9,7,0.5,N=100),
getprob(9,8,0.5,N=100),
getprob(9,9,0.5,N=100)
)

rbind(
getprob(9,1,0.2,N=100),
getprob(9,2,0.2,N=100),
getprob(9,3,0.2,N=100),
getprob(9,4,0.2,N=100),
getprob(9,5,0.2,N=100),
getprob(9,6,0.2,N=100),
getprob(9,7,0.2,N=100),
getprob(9,8,0.2,N=100),
getprob(9,9,0.2,N=100)
)
rbind(
getprob(9,1,0.9,N=100),
getprob(9,2,0.9,N=100),
getprob(9,3,0.9,N=100),
getprob(9,4,0.9,N=100),
getprob(9,5,0.9,N=100),
getprob(9,6,0.9,N=100),
getprob(9,7,0.9,N=100),
getprob(9,8,0.9,N=100),
getprob(9,9,0.9,N=100)
)

r_try <- 1:9
x_try <- seq(from = 0.1, to = 1, by = 0.1)
jj <- as.matrix(expand.grid(r_try, x_try))
II <- t(apply(jj, 1, function(v){getprob(9, v[1], v[2], N=1e4)}))
head(II)

r_true <- 3
x_true <- 0.85
o <- getprob(9, r_true, x_true, 17)  # observations
o[1:2]
o[-(1:2)]

support_func <- function(o, v){
    v <- v[-(1:2)]
    v <- v/sum(v)
    dmultinom(o[-(1:2)], prob=v, log=TRUE)
}
support <- rep(0,nrow(II))
for(i in seq_len(nrow(II))){
    support[i] <- support_func(o,II[i,])
}

support <- support-max(support)
support <- pmax(support, -10)
support  <- matrix(support,length(r_try),length(x_try))
contour(r_try, x_try, support, xlab="r", ylab="x", levels = -(1:5)*2)
abline(v = r_true)
abline(h = x_true)
abline(v = 1:9, col = 'gray', lty = 3, lwd = 0.3)
filled.contour(r_try, x_try, support, xlab = "r", ylab = "x")

```{R showjj} jj <- cbind(II[,1:2],c(support)) jj[jj[,2] == 1,]

#  Analysis from a single observation


We translate $\alpha\in [0,1]$ to a number from 1 to $n$ inclusive
with $\lceil n\alpha\rceil$

Write a likelihood function for 9,12 [that is, came 9th in a class of
12] There are 8 better and 3 worse than the focal student.


```r
single_like_fun <- function(n, a,          # n competitors, places a [e.g a=3 -> third]
                            a_try, x_try,  # values to try
                            N=1e2){        # N is number of trials
    p_try <- ceiling(n * a_try)
    jj <- as.matrix(expand.grid(p_try, x_try))
    II <- t(apply(jj,1,function(v){getprob(n, v[1], v[2], N)}))
    matrix(II[, ceiling(a*n) + 2], length(a_try), length(x_try))
}

a_try <- seq(from = 0.01, to = 0.99, len = 10)
x_try <- seq(from = 0.2, to = 1, by = 0.1)
allsupp <- list(
    log(single_like_fun(12, a =  8/12, a_try, x_try)),  # pure
    log(single_like_fun(17, a = 12/17, a_try, x_try)),  # pure
    log(single_like_fun(23, a = 18/23, a_try, x_try)),  # pure
    log(single_like_fun(14, a =  5/14, a_try, x_try)),  # applied
    log(single_like_fun(13, a =  3/13, a_try, x_try)),  # applied
    log(single_like_fun(11, a =  2/11, a_try, x_try)),  # applied
    log(single_like_fun(13, a =  3/13, a_try, x_try))   # applied
)

S <- Reduce("+",allsupp)
S <- pmax(S-max(S),-10)
contour(a_try, x_try, S, levels = -1*(1:7))
jj <- which(S == max(S), arr.ind = TRUE)
points(a_try[jj[1]], x_try[jj[2]], pch=16, col = 'red')

S_pure    <- Reduce("+", allsupp[1:3])
S_applied <- Reduce("+", allsupp[4:7])
S_all     <- Reduce("+", allsupp[1:7])
max(S_all)
max(S_pure)
max(S_applied)
max(S_pure) + max(S_applied)
max(S_pure) + max(S_applied) - max(S_all)
pchisq(max(S_pure) + max(S_applied) - max(S_all),df=2,lower.tail=FALSE)

Parkrun

In the following chunks, each execution of single_like_fun(...,N=1000) takes about 5 minutes

howmany <- 1000

print(system.time(parkrun01 <- single_like_fun(259, 177/259, N = howmany, a_try, x_try)))

parkrun01

 print(system.time(parkrun02 <- single_like_fun(305, 222/305, N = howmany, a_try, x_try)))

parkrun02

print(system.time(parkrun03 <- single_like_fun(297, 206/297, N = howmany, a_try, x_try)))

parkrun03

print(system.time(parkrun04 <- single_like_fun(241, 142/241, N = howmany, a_try, x_try)))

parkrun04

print(system.time(parkrun05 <- single_like_fun(179, 118/179, N = howmany, a_try, x_try)))

parkrun05

print(system.time(parkrun06 <- single_like_fun(338, 224/338, N = howmany, a_try, x_try)))

parkrun06

print(system.time(parkrun07 <- single_like_fun(203, 128/203, N = howmany, a_try, x_try)))

parkrun07

c(
sum(parkrun01),
sum(parkrun02),
sum(parkrun03),
sum(parkrun04),
sum(parkrun05),
sum(parkrun06),
sum(parkrun07)
)

contour(a_try,x_try,parkrun05)
contour(a_try,x_try,parkrun06)
parkrun01*parkrun02*parkrun03*parkrun04*parkrun05*parkrun06*parkrun07

Formula 1

n <- 20 # 20 drivers
multiple_like_fun <- function(n,           # n competitors
                            a_try, x_try,  # values to try
                            N=1e2){        # N is number of trials
    p_try <- ceiling(n * a_try)
    jj <- as.matrix(expand.grid(p_try, x_try))
    t(apply(jj,1,function(v){getprob(n, v[1], v[2], N)}))
}

o <- c(6,8,3,5,2,1,2,4,2,2,20,3,2,5,1,3,4,1,4,2,6,6,10) # lando norris 2024

II <- multiple_like_fun(20, a_try,x_try,1e4)
head(II)
nrow(II)
length(a_try)
length(x_try)

like <- list()
dim(II)
o
for(i in seq_along(o)){
    like[[i]] <-  matrix(II[, o[i] + 2], length(a_try), length(x_try))
}

like
alllike <- Reduce("*",like)
alllike 
alllike/max(alllike)
S <- log(alllike)
S <- S-max(S)
S <- pmax(S,-15)
contour(S,levels = -(0:6)*2, lwd=c(1,4,1,1,1,1,1))
filled.contour(S)