International Mathematical Olympiad using Plackett-Luce
robin hankin
2025-04-07
R Markdown
oecd <- as.character(read.csv("oecd.txt",header=FALSE)$V1)
oecd## [1] "Australia"
## [2] "Austria"
## [3] "Belgium"
## [4] "Canada"
## [5] "Chile"
## [6] "Colombia"
## [7] "Costa Rica"
## [8] "Czech Republic"
## [9] "Denmark"
## [10] "Estonia"
## [11] "Finland"
## [12] "France"
## [13] "Germany"
## [14] "Greece"
## [15] "Hungary"
## [16] "Iceland"
## [17] "Ireland"
## [18] "Israel"
## [19] "Italy"
## [20] "Japan"
## [21] "Democratic People's Republic of Korea"
## [22] "Latvia"
## [23] "Lithuania"
## [24] "Luxembourg"
## [25] "Mexico"
## [26] "Netherlands"
## [27] "New Zealand"
## [28] "Norway"
## [29] "Poland"
## [30] "Portugal"
## [31] "Slovakia"
## [32] "Slovenia"
## [33] "Spain"
## [34] "Sweden"
## [35] "Switzerland"
## [36] "Turkiye"
## [37] "United Kingdom"
## [38] "United States of America"oecd <- oecd[!(oecd %in% c("Chile", "Democratic People's Republic of Korea", "Turkiye"))]
years <- 2001:2024
jj <- read.table("imo.txt", header=TRUE)
ranks <- list()
for(i in seq_along(years)){
ranks[[i]] <- subset(jj, jj$year == years[[i]])$country
}
names(ranks) <- yearsfor(i in seq_along(ranks)){
jj <- ranks[[i]]
ranks[[i]] <- jj[jj %in% oecd]
}uk <- seq_along(ranks)*0
for(i in seq_along(ranks)){
uk[i] <- which(ranks[[i]] == "Australia")
}
uk## [1] 8 10 10 10 10 12 7 5 8 7 10 8 6 4 2 10 13 3 6 4 9 11 10 16BT <- 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
}
getprob <- function(n, r, x, N=1e2){
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)))
}
n <- length(oecd)
a_try <- seq(from = 0.01, to = 0.99, len = 20)
x_try <- seq(from = 0.2, to = 1, by = 0.05)
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)}))
}print(system.time(II <- multiple_like_fun(length(oecd), a_try, x_try, 1e4)))## user system elapsed
## 909.913 0.025 910.026head(II)## r.Var1 x.Var2 p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12
## [1,] 1 0.2 7977 1866 152 5 0 0 0 0 0 0 0 0
## [2,] 3 0.2 333 1339 6318 1871 136 3 0 0 0 0 0 0
## [3,] 4 0.2 76 257 1349 6414 1767 133 4 0 0 0 0 0
## [4,] 6 0.2 3 14 48 283 1283 6416 1802 147 4 0 0 0
## [5,] 8 0.2 0 0 4 7 62 241 1335 6470 1729 151 1 0
## [6,] 10 0.2 0 0 0 0 2 4 63 270 1267 6434 1808 150
## p13 p14 p15 p16 p17 p18 p19 p20 p21 p22 p23 p24 p25 p26 p27 p28 p29 p30
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## p31 p32 p33 p34 p35
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0
## [6,] 0 0 0 0 0nrow(II)## [1] 340length(a_try)## [1] 20length(x_try)## [1] 17like <- list()
dim(II)## [1] 340 37o <- uk # o for observation
o## [1] 8 10 10 10 10 12 7 5 8 7 10 8 6 4 2 10 13 3 6 4 9 11 10 16for(i in seq_along(o)){
like[[i]] <- matrix(II[, o[i] + 2], length(a_try), length(x_try))
}alllike <- Reduce("*",like)
S <- log(alllike)
S <- S-max(S)
S <- pmax(S, -9.999)
jj <- S
jj[jj < -9.99] <- NA
round(jj,2)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA NA NA NA NA NA NA -6.59
## [4,] NA NA NA NA NA NA NA NA NA NA NA -3.15 -2.05
## [5,] NA NA NA NA NA NA NA NA NA -1.63 -0.51 0.00 -1.80
## [6,] NA NA NA NA NA NA NA -9.24 -6.62 -4.37 -4.08 -4.36 -5.99
## [7,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [8,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [9,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [10,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [11,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [12,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [13,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [14,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [15,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [16,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [17,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [18,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [19,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [20,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [,14] [,15] [,16] [,17]
## [1,] NA -6.33 NA NA
## [2,] -5.79 -6.63 NA NA
## [3,] -4.58 -6.34 NA NA
## [4,] -3.78 -7.66 NA NA
## [5,] -5.32 -9.38 NA NA
## [6,] -9.10 NA NA NA
## [7,] NA NA NA NA
## [8,] NA NA NA NA
## [9,] NA NA NA NA
## [10,] NA NA NA NA
## [11,] NA NA NA NA
## [12,] NA NA NA NA
## [13,] NA NA NA NA
## [14,] NA NA NA NA
## [15,] NA NA NA NA
## [16,] NA NA NA NA
## [17,] NA NA NA NA
## [18,] NA NA NA NA
## [19,] NA NA NA NA
## [20,] NA NA NA NAcontour(a_try, x_try, S,
xlab = "BT strength (scaled)", ylab = "x",
levels = -(0:6)*2, lwd=c(1,4,1,1,1,1,1))
dim(S)## [1] 20 17contour(a_try[1:8], x_try[6:17], S[1:8, 6:17], levels = -(0:6)*2, lwd=c(1,4,1,1,1,1,1))
filled.contour(S)
like[[1]]## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 0 0 0 0 0 0 0 0 3 19 66 212 445
## [2,] 0 0 0 0 2 3 11 41 105 217 354 605 680
## [3,] 0 0 2 2 19 40 106 211 344 526 695 779 820
## [4,] 147 266 442 653 902 1056 1260 1429 1437 1377 1244 1059 890
## [5,] 6470 5535 4872 4241 3652 3118 2659 2167 1859 1509 1278 986 822
## [6,] 270 404 495 623 702 794 817 822 856 857 743 681 614
## [7,] 8 23 50 83 121 174 226 319 334 376 399 449 454
## [8,] 3 7 20 30 29 81 107 148 186 266 313 347 361
## [9,] 0 0 3 5 4 13 26 37 71 112 163 188 245
## [10,] 0 0 0 1 0 1 7 13 32 51 75 114 158
## [11,] 0 0 0 0 0 0 1 3 9 22 39 63 99
## [12,] 0 0 0 0 0 0 0 2 4 17 12 30 75
## [13,] 0 0 0 0 0 0 0 0 0 4 13 15 53
## [14,] 0 0 0 0 0 0 0 1 1 1 7 17 32
## [15,] 0 0 0 0 0 0 0 0 0 1 2 1 28
## [16,] 0 0 0 0 0 0 0 0 0 0 2 3 16
## [17,] 0 0 0 0 0 0 0 0 0 2 2 4 12
## [18,] 0 0 0 0 0 0 0 0 0 0 1 2 5
## [19,] 0 0 0 0 0 0 0 0 0 0 1 2 5
## [20,] 0 0 0 0 0 0 0 0 0 0 0 0 9
## [,14] [,15] [,16] [,17]
## [1,] 617 593 434 277
## [2,] 643 606 446 233
## [3,] 737 583 442 288
## [4,] 716 562 409 260
## [5,] 658 476 432 282
## [6,] 517 462 365 320
## [7,] 404 416 345 301
## [8,] 341 379 343 275
## [9,] 267 300 279 298
## [10,] 212 263 305 271
## [11,] 148 227 274 261
## [12,] 108 181 247 265
## [13,] 80 133 223 254
## [14,] 57 145 173 256
## [15,] 42 94 185 266
## [16,] 31 85 191 295
## [17,] 33 69 180 253
## [18,] 17 61 135 276
## [19,] 13 64 150 299
## [20,] 15 49 120 298Some thoughts about including variable numbers of competitors
Above we used OECD countries for Australia to compete against. Here I use all the competitors:
a <- read.table("imo.txt", header = TRUE)
min(a$year)## [1] 2001max(a$year)## [1] 2024which(subset(a,a$year == 2015)$country == "Australia")## [1] 6competitors <- function(year,a){ table(a$year)[as.character(year)] }
competitors(2013, a)## 2013
## 97rank <- function(place, y, a){
which(subset(a,a$year == y)$country == place)
}
rank("Australia",2001,a)## [1] 25rank("Australia",2002,a)## [1] 26rank("Australia",2024,a)## [1] 39cbind(place = sapply(2001:2024,function(y){rank("Australia",y,a)}),
competitors = table(a$year)
)## place competitors
## 2001 25 83
## 2002 26 84
## 2003 26 82
## 2004 27 85
## 2005 25 91
## 2006 26 90
## 2007 22 93
## 2008 19 97
## 2009 23 104
## 2010 15 96
## 2011 25 101
## 2012 27 100
## 2013 15 97
## 2014 11 101
## 2015 6 104
## 2016 25 109
## 2017 34 111
## 2018 11 107
## 2019 18 112
## 2020 8 105
## 2021 18 107
## 2022 29 104
## 2023 23 112
## 2024 39 108Now try New Zealand:
cbind(place = sapply(2001:2024,function(y){rank("New Zealand",y,a)}),
competitors = table(a$year)
)## place competitors
## 2001 44 83
## 2002 36 84
## 2003 57 82
## 2004 58 85
## 2005 38 91
## 2006 58 90
## 2007 51 93
## 2008 72 97
## 2009 66 104
## 2010 29 96
## 2011 30 101
## 2012 53 100
## 2013 50 97
## 2014 60 101
## 2015 50 104
## 2016 53 109
## 2017 47 111
## 2018 46 107
## 2019 58 112
## 2020 47 105
## 2021 70 107
## 2022 58 104
## 2023 66 112
## 2024 32 108