exponential Bradley-Terry: likelihood-based methods on a grid of values
In hyper2: The Hyperdirichlet Distribution, Mark 2

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

(takes maybe twenty minutes to run without cach at $N=9$)

First we define N, the number of competitors. This is used consistently in the whole Rmd file, except for the first few bits where I hard-code five or 11 competitors for ease of exposition. In principle, the Rmd file should process whatever the value of N.

N <- 9  # number of competitors

Now we define BT() which specifies Bradley-Terry strengths $\alpha x^i$:

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) <- str_pad(1:n, width=ceiling(log10(n)), pad = "0")
    out
}
BT(4,0.5)
sum(BT(4,0.5))

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

makeracetable <- function(no_of_races, no_of_competitors, x){
    t(replicate(no_of_races,
                as.numeric(rrace(BT(no_of_competitors, x)))))
}

makeracetable(no_of_races=10, no_of_competitors=5, x=0.7)

get_one_competitor <- function(racetable, competitor){
    jj <- apply(racetable, 1, function(x){which(x==competitor)})
    tabulate(jj, ncol(racetable))
}

X <- makeracetable(1e4,11,0.7)  # using N=11 here
get_one_competitor(X,4)
get_one_competitor(X,11)

rbind(
get_one_competitor(X,1) ,
get_one_competitor(X,2) ,
get_one_competitor(X,9) ,
get_one_competitor(X,10)
)

t(sapply(seq_len(11),function(i){get_one_competitor(X,i)}))

OK now to put above code in to functional form:

get_all_competitors <- function(no_of_races, no_of_competitors, x){
    X <- makeracetable(no_of_races=no_of_races, no_of_competitors=no_of_competitors, x=x)
    out <- t(sapply(seq_len(no_of_competitors),function(i){get_one_competitor(X,i)}))
    dimnames(out) <- list(
        true_rank = paste("r_", seq_len(no_of_competitors), sep=""),
        observed_rank = paste("c_", seq_len(no_of_competitors), sep=""))
    return(out/no_of_races)
}

Now use it:

get_all_competitors(100, N, 0.5)

Above, we see from the first line that, with a sample of 100 races, the top competitor comes first with (estimated) probability about 0.XXXX, second with probability abot 0.XXX, and so on, coming fifth with probability 0.XXXXX

get_all_competitors(100, N, 0.99)

Above, all competitors are about the same. We see a more uniform spread of observations.

get_all_competitors(100, N, 0.1)

Above, each competitor is ten times as good as the next one. We see underdispersed results. So what we need to do now is a "meaningful grid". This is somewhat computationally intensive.

no_of_races <- 1000
no_of_competitors <- N

true_x <- seq(from=0.1, by=0.1, to=0.9)

First, we combine a whole bunch of X matrices together:

allX <- function(no_of_competitors, no_of_races, true_x){
    out <- NULL
    for(i in seq_along(true_x)){
        out <- rbind(out,
                     cbind(
                         true_x[i],
                         seq_len(no_of_competitors),
                         get_all_competitors(
             no_of_races = no_of_races,
             no_of_competitors = no_of_competitors,
             true_x[i])
                     )
                     )
    }
    rownames(out) <- NULL
    jj <- colnames(out)
    jj[1:2] <- c("x","r")
    colnames(out) <- jj
    return(out)
}

allX9 <- allX(
          no_of_competitors = N,
          no_of_races = 50000,
          true_x = seq(from=0.1, to=0.9, by=0.1)
      )

head(allX9)
allX9[seq(from=round(nrow(allX9)/2 - 5), to=round(nrow(allX9)/2 + 5)),]
tail(allX9)

Now some observations:

make_obs <- function(no_of_observations, no_of_competitors, r_true, x_true){
    kk <- makeracetable(no_of_races = no_of_observations, no_of_competitors = no_of_competitors, x = x_true)
    obs <- tabulate(apply(kk,1,function(o){which(o == r_true)}), nbins=no_of_competitors)
    return(obs)
}

make_obs(no_of_observations=20, no_of_competitors=N, x_true = 0.7, r_true = 6)
make_obs(no_of_observations=20, no_of_competitors=N, x_true = 0.7, r_true = 9)
make_obs(no_of_observations=20, no_of_competitors=N, x_true = 0.1, r_true = 6)


makemaxliketable <- function(n, nobs, no_of_competitors, allX, r_true, x_true){

    maxlike <- function(obs){
        support <- function(i){dmultinom(obs, prob=allX[i, -(1:2)], log=TRUE)}
        p <- sapply(seq_len(nrow(allX)), support)
        allX[which.max(p), 1:2]
    }

    evaluate <- matrix(NA,n,2)
    for(i in seq_len(n)){
        obs <- make_obs(no_of_observations=nobs, no_of_competitors=no_of_competitors, x_true=x_true, r_true=r_true)
        evaluate[i,] <- maxlike(obs)
    }
    return(table(x=evaluate[,1], r=evaluate[,2]))
}

biasx <- function(liketable, x_true){
    values <- as.numeric(rownames(liketable))
    probs <- rowSums(liketable) / sum(liketable) # sum(probs) == 1
    sum(probs*(values - x_true))
}

biasr <- function(liketable, r_true){
    values <- as.numeric(colnames(liketable))
    probs <- colSums(liketable) / sum(liketable) # sum(probs) == 1
    sum(probs*(values - r_true))

}

msex <- function(liketable, x_true){
    values <- as.numeric(rownames(liketable))
    probs <- rowSums(liketable) / sum(liketable) # sum(probs) == 1
    sum(probs*(values - x_true)^2)
}

mser <- function(liketable, r_true){
    values <- as.numeric(colnames(liketable))
    probs <- colSums(liketable) / sum(liketable) # sum(probs) == 1
    sum(probs*(values - r_true)^2)
}

(jj <- makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.8))
biasx(jj, x_true = 0.8)
biasr(jj, r_true = 4)
msex(jj, x_true = 0.8)
mser(jj, r_true = 4)

Vary nothing

Just do the same thing a few times to gauge the variability in the system:

jj1 <- makemaxliketable(n=100, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj2 <- makemaxliketable(n=100, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj3 <- makemaxliketable(n=100, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj4 <- makemaxliketable(n=100, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj5 <- makemaxliketable(n=100, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj6 <- makemaxliketable(n=100, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj7 <- makemaxliketable(n=100, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj8 <- makemaxliketable(n=100, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj9 <- makemaxliketable(n=100, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj1
jj2

```r
fish <-
rbind(
c(biasr(jj1, 4), mser(jj1, 3.8), biasx(jj1, 0.7), msex(jj1, 0.7)),
c(biasr(jj2, 4), mser(jj2, 3.8), biasx(jj2, 0.7), msex(jj2, 0.7)),
c(biasr(jj3, 4), mser(jj3, 3.8), biasx(jj3, 0.7), msex(jj3, 0.7)),
c(biasr(jj4, 4), mser(jj4, 3.8), biasx(jj4, 0.7), msex(jj4, 0.7)),
c(biasr(jj5, 4), mser(jj5, 3.8), biasx(jj5, 0.7), msex(jj5, 0.7)),
c(biasr(jj6, 4), mser(jj6, 3.8), biasx(jj6, 0.7), msex(jj6, 0.7)),
c(biasr(jj7, 4), mser(jj7, 3.8), biasx(jj7, 0.7), msex(jj7, 0.7)),
c(biasr(jj8, 4), mser(jj8, 3.8), biasx(jj8, 0.7), msex(jj8, 0.7)),
c(biasr(jj9, 4), mser(jj9, 3.8), biasx(jj9, 0.7), msex(jj9, 0.7))
)
colnames(fish) <- c("bias(r)", "mse(r)", "bias(x)", "mse(x)")

fish
mean(fish[,1])
mean(fish[,2])
sd(fish[,1])
mean(fish[,2])/sqrt(100)

jj1 <- makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj2 <- makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj3 <- makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj4 <- makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj5 <- makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj6 <- makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj7 <- makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj8 <- makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj9 <- makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
jj1
jj2

```r
fish <-
rbind(
c(biasr(jj1, 4), mser(jj1, 3.8), biasx(jj1, 0.7), msex(jj1, 0.7)),
c(biasr(jj2, 4), mser(jj2, 3.8), biasx(jj2, 0.7), msex(jj2, 0.7)),
c(biasr(jj3, 4), mser(jj3, 3.8), biasx(jj3, 0.7), msex(jj3, 0.7)),
c(biasr(jj4, 4), mser(jj4, 3.8), biasx(jj4, 0.7), msex(jj4, 0.7)),
c(biasr(jj5, 4), mser(jj5, 3.8), biasx(jj5, 0.7), msex(jj5, 0.7)),
c(biasr(jj6, 4), mser(jj6, 3.8), biasx(jj6, 0.7), msex(jj6, 0.7)),
c(biasr(jj7, 4), mser(jj7, 3.8), biasx(jj7, 0.7), msex(jj7, 0.7)),
c(biasr(jj8, 4), mser(jj8, 3.8), biasx(jj8, 0.7), msex(jj8, 0.7)),
c(biasr(jj9, 4), mser(jj9, 3.8), biasx(jj9, 0.7), msex(jj9, 0.7))
)
colnames(fish) <- c("bias(r)", "mse(r)", "bias(x)", "mse(x)")

fish
mean(fish[,1])
mean(fish[,2])
sd(fish[,1])
mean(fish[,2])/sqrt(1000)

Above we have that the bias of $r$ has a standard deviation, given (approximately at least) by the standard error of the mean.

Vary `x_true`

vxt <- list(
makemaxliketable(n=4000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.1),
makemaxliketable(n=4000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.2),
makemaxliketable(n=4000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.3),
makemaxliketable(n=4000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.4),
makemaxliketable(n=4000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.5),
makemaxliketable(n=4000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.6),
makemaxliketable(n=4000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7),
makemaxliketable(n=4000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.8),
makemaxliketable(n=4000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.9))

jj <- rbind(
c(biasr(vxt[[1]], 4), biasx(vxt[[1]], 0.1), mser(vxt[[1]], 4), msex(vxt[[1]], 0.1)),
c(biasr(vxt[[2]], 4), biasx(vxt[[2]], 0.2), mser(vxt[[2]], 4), msex(vxt[[2]], 0.2)),
c(biasr(vxt[[3]], 4), biasx(vxt[[3]], 0.3), mser(vxt[[3]], 4), msex(vxt[[3]], 0.3)),
c(biasr(vxt[[4]], 4), biasx(vxt[[4]], 0.4), mser(vxt[[4]], 4), msex(vxt[[4]], 0.4)),
c(biasr(vxt[[5]], 4), biasx(vxt[[5]], 0.5), mser(vxt[[5]], 4), msex(vxt[[5]], 0.5)),
c(biasr(vxt[[6]], 4), biasx(vxt[[6]], 0.6), mser(vxt[[6]], 4), msex(vxt[[6]], 0.6)),
c(biasr(vxt[[7]], 4), biasx(vxt[[7]], 0.7), mser(vxt[[7]], 4), msex(vxt[[7]], 0.7)),
c(biasr(vxt[[8]], 4), biasx(vxt[[8]], 0.8), mser(vxt[[8]], 4), msex(vxt[[8]], 0.8)),
c(biasr(vxt[[9]], 4), biasx(vxt[[9]], 0.9), mser(vxt[[9]], 4), msex(vxt[[9]], 0.9))
)
colnames(jj) <- c("bias(r)","bias(x)", "mse(r)", "mse(x)")
rownames(jj) <- seq(from=0.1,to=0.9,by=0.1)
jj
matplot(jj,type="b")
cbind(`bias(r)`=jj[,1], `SE(bias)` = jj[,3]/sqrt(4000))
cbind(`bias(x)`=jj[,2], `SE(bias)` = jj[,4]/sqrt(4000))

Above we see that varying x_true from 0.1 to 0.9 increases the variance of $r_\text{estimated}$

Vary `r_true`

VR <- list(
makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 1, x_true = 0.7),
makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 2, x_true = 0.7),
makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 3, x_true = 0.7),
makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7),
makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 5, x_true = 0.7),
makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 6, x_true = 0.7),
makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 7, x_true = 0.7),
makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 8, x_true = 0.7),
makemaxliketable(n=1000, nobs=10, no_of_competitors=N, allX=allX9, r_true = 9, x_true = 0.7)
)

jj <- rbind(
c(biasr(VR[[1]], 1), biasx(VR[[1]], 0.7), mser(VR[[1]], 1), msex(VR[[1]], 0.7)),
c(biasr(VR[[2]], 2), biasx(VR[[2]], 0.7), mser(VR[[2]], 2), msex(VR[[2]], 0.7)),
c(biasr(VR[[3]], 3), biasx(VR[[3]], 0.7), mser(VR[[3]], 3), msex(VR[[3]], 0.7)),
c(biasr(VR[[4]], 4), biasx(VR[[4]], 0.7), mser(VR[[4]], 4), msex(VR[[4]], 0.7)),
c(biasr(VR[[5]], 5), biasx(VR[[5]], 0.7), mser(VR[[5]], 5), msex(VR[[5]], 0.7)),
c(biasr(VR[[6]], 6), biasx(VR[[6]], 0.7), mser(VR[[6]], 6), msex(VR[[6]], 0.7)),
c(biasr(VR[[7]], 7), biasx(VR[[7]], 0.7), mser(VR[[7]], 7), msex(VR[[7]], 0.7)),
c(biasr(VR[[8]], 8), biasx(VR[[8]], 0.7), mser(VR[[8]], 8), msex(VR[[8]], 0.7)),
c(biasr(VR[[9]], 9), biasx(VR[[9]], 0.7), mser(VR[[9]], 9), msex(VR[[9]], 0.7))
)
colnames(jj) <- c("bias(r)","bias(x)", "mse(r)", "mse(x)")
rownames(jj) <- seq(from=0.1,  to=0.9, by=0.1)
jj
matplot(jj,type="b")
cbind(`bias(r)`=jj[,1], `SE(bias)` = jj[,3]/sqrt(1000))
cbind(`bias(x)`=jj[,2], `SE(bias)` = jj[,4]/sqrt(1000))

Vary `nobs`

VN <- list(
makemaxliketable(n=1000, nobs= 10, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7),
makemaxliketable(n=1000, nobs= 20, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7),
makemaxliketable(n=1000, nobs= 50, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7),
makemaxliketable(n=1000, nobs=100, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7),
makemaxliketable(n=1000, nobs=200, no_of_competitors=N, allX=allX9, r_true = 4, x_true = 0.7)
)

jj <- rbind(
c(biasr(VN[[1]], 4), biasx(VN[[1]], 0.7), mser(VN[[1]], 4), msex(VN[[1]], 0.7)),
c(biasr(VN[[2]], 4), biasx(VN[[2]], 0.7), mser(VN[[2]], 4), msex(VN[[2]], 0.7)),
c(biasr(VN[[3]], 4), biasx(VN[[3]], 0.7), mser(VN[[3]], 4), msex(VN[[3]], 0.7)),
c(biasr(VN[[4]], 4), biasx(VN[[4]], 0.7), mser(VN[[4]], 4), msex(VN[[4]], 0.7)),
c(biasr(VN[[5]], 4), biasx(VN[[5]], 0.7), mser(VN[[5]], 4), msex(VN[[5]], 0.7))
)
colnames(jj) <- c("bias(r)","bias(x)", "mse(r)", "mse(x)")
rownames(jj) <- c(10,20,50,100,200)
jj
matplot(jj,type="b")
cbind(`bias(r)`=jj[,1], `SE(bias)` = jj[,3]/sqrt(1000))
cbind(`bias(x)`=jj[,2], `SE(bias)` = jj[,4]/sqrt(1000))