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)))
}

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.

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)

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

howmany <- 100

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*10, 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

M <- matrix(c(
259, 177,
305, 222,
297, 206,
241, 142,
179, 118,
338, 224,
203, 128),ncol=2,byrow=TRUE)
M
plot(M[,2]/M[,1])

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