In RobinHankin/hyper2: The Hyperdirichlet Distribution, Mark 2

set.seed(1)
knitr::opts_chunk$set(echo = TRUE)
library("hyper2")

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

To cite the hyper2 package in publications, please use @hankin2017_rmd. Here I analyse results from Curacao 1962 using hyper3 formalism using the weighted Draw monster approach outlined in kka.Rmd:

$$ \left(\frac{\lambda p_1, D(p_1+p_2),p_2}{\lambda p_1+D(p_1+p_2)+p_2}\right)^{(a,b,c)} \left(\frac{p_1, D(p_1+p_2),\lambda p_2}{p_1+D(p_1+p_2)+\lambda p_2}\right)^{(d,e,f)} $$

\newcommand{\wld}[3]{\mathord{+}#1\,\mathord{=}#2\,\mathord{-}#3}

Above, $p_1$ plays $p_2$ a total of $a+b+c+d+e+f$ times; we have $a+b+c$ matches with $p_1$ playing white and $p_2$ black, with $a$ wins, $b$ draws and $c$ losses [write $\wld{a}{b}{c}$], and $d+e+f$ matches with $p_1$ playing black and $p_2$ playing white with $\wld{d}{e}{f}$.

Short story

With this model [discussed in a broader context at kka.Rmd], whatever the value of $D$ and $\lambda$, there is a pathology which BT cannot deal with. Still not quite 100% sure but it seems that the maximum likelihood estimate of the players' strengths is zero except for Fisher and Petrosian; and indeed we see Fisher at 57% and Petrosian at 42% [with $\lambda = 1.1, D=1.888$]. This on the grounds that these players never lost playing black (there are some $2\times 3$ tables at the end). This might be considered a defect of the probability model.

Longer story

cur_matches <- read.table("stockholm1962_matches.txt",header=FALSE)
colnames(cur_matches) <- c("white","black","result")
head(cur_matches)
nrow(cur_matches)

Above, we see on the first line that Aaron (white) lost to Barcza (black). We will set up hyper2 and hyper3 objects corresponding to likelihood functions for these observations.

jj <- read.table("stockholm1962.txt",header=FALSE)[,1:2]
(players <- jj$V1)
(nationalities <- jj$V2)

f <- function(lambda,D){
  H3 <- hyper3()
  for(i in seq_len(nrow(cur_matches))){
    white_player <- cur_matches[i,1]
    black_player <- cur_matches[i,2]
    result <- cur_matches[i,3]
    if(result == "1-0"){ # white player wins
       num <- lambda 
       names(num) <- white_player
    } else if(result == "0-1"){  # black player wins
       num <- 1
       names(num) <- black_player
    } else if (result == "1/2-1/2"){ # draw, duh
      num <- c(D,D)
      names(num) <- c(white_player,black_player)
    } else {
      stop("this cannot happen")
    }

    H3[num] %<>% inc()       
    den <- c(lambda+D,1+D)
    names(den) <- c(white_player,black_player)
    H3[den] %<>% dec()
  }
  return(H3)
}

Now use it

H3 <- f(1.1,1.888)

H3
summary(H3)

equalp.test(H3,n=1)

m1121 <- maxp(f(1.1,2.1),n=1,give=TRUE)
m1123 <- maxp(f(1.1,2.3),n=1,give=TRUE)

pie(m1121$par)

A break here.

profl <- function(v){
    lambda <- v[1]
    D <- v[2]
    H <- f(lambda,v)
    return(loglik(maxp(H,n=1),H))
}

Now apply profl():

profl(c(1.1,2.1))
#profl(c(1.2,2.1))

np <- 4
lambda <- seq(from=0.1,to=1.5,len=np)
D <- seq(from=0.1,to=0.4,len=np)
jj <- as.matrix(expand.grid(lambda,D))
LL <- apply(jj,1,profl)
LL

LL <- LL - max(LL)
m <- matrix(LL,np,np)
m
filled.contour(lambda,D,m,xlab="lambda",ylab="D",nlevels=40)

res <- cur_matches$V3
res
table(res)
white_wins <- cur_matches[res=="1-0"    ,]
black_wins <- cur_matches[res=="0-1"    ,]
draws      <- cur_matches[res=="1/2-1/2",]
(players <- unique(sort(c(cur_matches[,1],cur_matches[,2]))))
jj <- rep(0,length(players))
names(jj) <- players
jj
plays_white_wins <- jj
plays_white_lose <- jj
plays_white_draw <- jj
plays_black_wins <- jj
plays_black_lose <- jj
plays_black_draw <- jj
for(i in seq_len(nrow(cur_matches))){
  white_player <- cur_matches[i,1]
  black_player <- cur_matches[i,2]
  result       <- cur_matches[i,3]  

  if(result == "1-0"){
      plays_white_wins[white_player] %<>% inc
      plays_black_lose[black_player] %<>% inc
  } else if(result == "0-1"){
      plays_white_lose[white_player] %<>% inc
      plays_black_wins[black_player] %<>% inc
  } else if(result == "1/2-1/2"){
      plays_white_draw[white_player] %<>% inc
      plays_black_draw[black_player] %<>% inc
  }
}

plays_white_wins
plays_white_lose
plays_black_wins
plays_black_lose

Now we need some way to describe each player's wins/losses playing white and black:

f <- function(player){
out <- matrix(c(
  plays_white_wins[player], plays_white_lose[player], plays_white_draw[player],
  plays_black_wins[player], plays_black_lose[player], plays_black_draw[player]
  ),byrow=TRUE,2,3)
  dimnames(out) <- list(colour=c("white","black"),match=c("win","lose","draw"))
  return(out)
}
f("Fischer")
f("Korchnoi")
f("Petrosian")
M <- f("Aaron")*0
for(p in players){ M <- M + f(p) }
M