In RobinHankin/hyper2: The Hyperdirichlet Distribution, Mark 2

knitr::opts_chunk$set(echo = TRUE)

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

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

@paul1995 state that "underdispersion will be seen very rarely in real-life situations".

library("hyper2")

f <- function(n,alpha){
  jj <- 1/alpha^(1:n)
  return(jj/sum(jj))
}

getlike <- function(v){
  alpha <- v[1]
  beta <- v[2]
 c(
tabulate(replicate(10000,which(rrace(setNames(f(26,alpha),letters))==letters[beta])),nbins=26)
,alpha,beta)
}

alpha_try <- seq(from=1,to=6,len=20)
beta_try <- 1:26
calculate_from_scratch <- FALSE
if(calculate_from_scratch){
    M <- as.matrix(expand.grid(alpha=alpha_try, beta=beta_try))
    system.time(L <- t(apply(M,1,getlike)))
    colnames(L) <- c(letters,'alpha','beta')
    saveRDS(L,file="likelihoodmatrix.Rdata")
} else {
    L <- readRDS("likelihoodmatrix.Rdata")
}
L_integer <- L  # counts
L[,1:26] <- L[,1:26]/mean(rowSums(L[,1:26]))
head(L_integer)
head(L)

colSums(L)
mL <- function(i){matrix(L[,i],length(alpha_try),length(beta_try))}
filled.contour(alpha_try,beta_try,mL(1))
filled.contour(alpha_try,beta_try,mL(5))
filled.contour(alpha_try,beta_try,mL(10))
filled.contour(alpha_try,beta_try,mL(15))
filled.contour(alpha_try,beta_try,mL(20))
filled.contour(alpha_try,beta_try,mL(26))

#plot(table(L_integer[,1:26]))

Now 10 in-silico races:

alphatrue <- 2.3
betatrue <- 11 
set.seed(1)
D <- replicate(17,which(rrace(setNames(f(26,alphatrue),letters))==letters[betatrue]))
sort(D)
plot(table(D))

likelihood <- apply(L[,D],1,prod) + 0.00
likelihood <- matrix(likelihood,length(alpha_try),length(beta_try))

round(alpha_try,2)
matplot(alpha_try,likelihood,type='b')
matplot(beta_try,t(likelihood),type='b')
matplot(alpha_try,log(likelihood)-log(max(likelihood)),type='b',ylim=c(-20,0))
matplot(beta_try,log(t(likelihood))-log(max(likelihood)),type='b',ylim=c(-20,0))
filled.contour(alpha_try,beta_try,(likelihood))
table(likelihood==0)
plot(f(26,alphatrue))
abline(v=betatrue)
contour(alpha_try,beta_try,likelihood)
abline(v=alphatrue)
abline(h=betatrue)
contour(alpha_try,beta_try,log(likelihood) - log(max(likelihood)),levels=-(0:10))
abline(v=alphatrue)
abline(h=betatrue)

\subsection{Park run}

"Parkrun" is a distributed community initiative that organises weekly timed 5 km runs/walks in parks worldwide~\cite{hindley2020}. A typical event will have 200 participants. Table~\ref{parkruntable} shows the author has completed a total of 15 parkruns to date, and from the first pair of numbers we see that 238 runners attended that particular parkrun, of whom the author placed 173.

\begin{table}[t] \caption{Parkrun results} \label{parkruntable} \begin{tabular}{ccccccccccc} \hline rank &173 & 165 & 172 & 199 & 181 & 229 & 177 & 222 & 206 & 142 \ runners&238 & 238 & 196 & 242 & 242 & 318 & 259 & 305 & 297 & 241 \ \ \ rank & 118 & 224 & 128 & 115 & 183 & & & & &\ runners& 179 & 338 & 203 & 245 & 254 & & & & &\ \hline \end{tabular} \end{table}

We may consider the author to have an unknown generalized Bradley Terry strength $a$: Figure~\ref{parkrunsupport} shows a support curve for $a$ with the evaluate, $\hat{a}\simeq 0.448$ shown. We also see two units of support~\cite{edwards1992} shown as a dotted line illustrating a support interval of about 0.315 - 0.567; we may be reasonably confident that the author's true BT strength lies in this range. In particular, note that a reasonable $H_0\colon a=\frac{1}{2}$ may not be rejected, the support for $H_0$ being only susceptible to very minor improvement [$\simeq 0.351$] by moving to the evaluate.

\begin{figure}[t]
\includegraphics[width=4in]{plotparkrun}  % plotparkrun.pdf made in very_simplified_likelihood.Rmd
\caption{Author's strength\label{parkrunsupport}}
\end{figure}

RobinHankin/hyper2 documentation built on April 21, 2024, 11:38 a.m.