In RobinHankin/hyper2: The Hyperdirichlet Distribution, Mark 2

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

Formula 1 grand prix racing, 2023 season.

setwd("~/rstudio/hyper2/inst")
finish <- read.table("formula1_2023.txt",header=TRUE)
finish <- finish[,seq_len(ncol(finish)-1)]
finish[,1:5]

grid <- read.table("formula1_starting_2023.txt",header=TRUE)
grid[,1:5]

replacenonzerowithorder <- function(x){
  nonzero <- x!=0
  x[nonzero] <- rank(x[nonzero])
  return(x)
}

order_non_NA_NZ <- function(x){ # leaves NA and 0 entries unchanged; replaces other elements with their order
  wanted <- !is.na(x)
  x[wanted] <- replacenonzerowithorder(x[wanted])
  return(x)
}

topn <- 15
finish <- finish[seq_len(topn),]
options(use_alabama=TRUE)

`makeH3` <- function(i,lambda,topn,finish,grid){
  suppressWarnings(a <- setNames(as.numeric(finish[,i]),rownames(finish)))
  suppressWarnings(b <- setNames(as.numeric(grid  [,i]),rownames(finish)))
  wanted <- !is.na(b)
  a <- a[wanted]
  b <- b[wanted]
  a[is.na(a)] <- 0 # did not finish
  a <- a[seq_len(topn)]
  a <- order_non_NA_NZ(a)
  H <- ordervec2supp(a)
  for(j in seq_along(b)){
    if(!is.na(names(b)[j])){H[names(b)[j]] <- as.weight(lambda^b[j])}
  }
  return(H)
}

likelam <- function(lambda,topn,finish,grid){
  H <- hyper3()
  for(i in seq_len(ncol(finish))){H <- H + makeH3(i,lambda=lambda,topn=topn,finish=finish,grid=grid)}
  return(H)
}

f <- function(lambda,topn,finish,grid){
  H <- likelam(lambda,topn,finish,grid)
  m <- maxp(H,n=1)
  return(c(loglik(m,H),lambda,m))
}

system.time(lam <- seq(from=0.85, to=1.01, len=15))
system.time(bigmat <- sapply(lam,f,topn=topn,finish=finish,grid=grid))

head(bigmat)

system.time(alpha_MLE <- optimize(function(a){f(a,topn=topn,finish=finish,grid=grid)[1]},c(0.90,0.95),maximum=TRUE))
alpha_MLE

like <- bigmat[1,]
like <- like - max(like)
plot(lam,like,type='b')
abline(h=c(0,-2))
abline(v=1)
segments(x0=0.9654,y0=-5,y1=0)
grid()

if(FALSE){print(alpha_MLE)}

H2 <- ordertable2supp(finish)
mH2 <- maxp(H2)

dotchart(mH2,pch=16)
dotchart(log(mH2),pch=16)

o <- likelam(lambda=0.9, topn=15,finish=finish,grid=grid)
loglik(equalp(o),o)

o <- likelam(lambda=0.95, topn=15,finish=finish,grid=grid)
loglik(equalp(o),o)

Above was 2023, now do the 2022 season

setwd("~/rstudio/hyper2/inst")
finish <- read.table("formula1_2022.txt",header=TRUE)
finish22 <- finish[,seq_len(ncol(finish)-1)]
finish22[,1:5]
grid22 <- read.table("formula1_starting_2022.txt",header=TRUE)
grid22[,1:5]

system.time(
    alpha_MLE22 <- optimize(function(a){f(a,topn=topn,finish=finish22,grid=grid22)[1]},c(0.95,0.99),maximum=TRUE)
    )
alpha_MLE22

lam <- seq(from=0.90, to=1.05, len=15)
system.time(bigmat <- sapply(lam,f,topn=topn,finish=finish22,grid=grid22))

head(bigmat)

like <- bigmat[1,]
like <- like - max(like)
plot(lam,like,type='b')
abline(h=c(0,-2))

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