knitr::opts_chunk$set(echo = TRUE) library("hyper2") library("pracma") library("magrittr")
knitr::include_graphics(system.file("help/figures/hyper2.png", package = "hyper2"))
Consider a race between a competitor of Bradley-Terry strength $a$ and $n$ competitors each of strength $b$; we require $a+b=1$. NB we have $n+1$ competitors in total, 1 of strength $a$ and $n$ of strength $b$. An observation is indexed by $r$, the number of $b$ clones finishing ahead of $a$, so $0\leqslant r\leqslant n$. The initial field strength is $a+nb=1+b(n-1)$.
Suppose our bro comes first (zero $b$'s ahead of him):
[ \mathcal{L}(a)=\frac{a}{a+nb}=\frac{1-b}{1+(n-1)b} ]
Comes second, one $b$ ahead of him:
[ \mathcal{L}(a)=\frac{b}{a+nb}\cdot\frac{a}{a+(n-1)b} ]
Third, two ahead of him:
[ \mathcal{L}(a)= \frac{b}{a+ n b}\cdot \frac{b}{a+(n-1)b}\cdot \frac{a}{a+(n-2)b} ]
Fourth, three ahead:
[ \mathcal{L}(a)= \frac{b}{a+ n b}\cdot \frac{b}{a+(n-1)b}\cdot \frac{b}{a+(n-2)b}\cdot \frac{a}{a+(n-3)b} ]
Placing after $r<n$ of the $n$ clones of strength $b$ finish:
[ \mathcal{L}(a)= \frac{b^ra}{ (a+nb)(a+(n-1)b)(a+(n-2)b)\cdots(a+(n-r)b)} = \frac{B-1}{(B+n-1)(B+n-2)\cdots(B+n-r-1)} ]
where $B=1/b$ [also using the fact that $a+b=1$].
This is
[ \mathcal{L}_{r,n}(a)=\frac{(B-1)(B+n-r-2)!}{(B+n-1)!},\qquad B=1/(1-a) ]
Then, using Stirling [viz $n!\sim n^ne^{-n}\sqrt{2\pi n}$] we get
[ \log(B-1) + ((B+n-r-2)\log(B+n-r-2)-(B+n-r-2)+\frac{\log(B+n-r-2)}{2}) - ((B+n -1)\log(B+n -1)-(B+n -1)+\frac{\log(B+n -1)}{2}) ]
If you give Sage this:
var('B n r') pretty_print_default(True) X=log(B-1) + ( + ((B+n-r-2)*log(B+n-r-2)-(B+n-r-2)+log(B+n-r-2)/2) - ((B+n -1)*log(B+n -1)-(B+n -1)+log(B+n -1)/2) ) X.taylor(n,Infinity,2)
It returns this:
[ -(r+1)\log n -\frac{(2B-3)r-r^2+2B-2}{2n} -\frac{3(2B-3)r^2-6B^2-6(B^2-3B+2)r+12B-5}{12n^2} +\log(B-1)+\mathcal{O}(n^{-3}) ]
Note that if we use instead $n~\sim n^ne^{-n}\sqrt{2\pi n}\left(1+\frac{1}{12n}\right)$, then only terms of $\mathcal{O}(n^{-2})$ and above are changed; the first-order terms and indeed the zeroth order terms are unaffected. Anyway, after simplifying and taking an asymptotic expansion to order $n^{-1}$, we find, to first order,
[ \mathcal{S}=\log\mathcal{L}_{r,n}(a)=\log(B-1)-B\frac{r+1}{n} + K + \mathcal{O}\left(n^{-2}\right) ]
(where $K$ is an arbitrary constant). Observing that $\partial^2\mathcal{S}/\partial B^2=-(B-1)^{-2}<0$, the evaluate would be unique; we have $\hat{B}=\frac{n+r+1}{r+1}+\mathcal{O}\left(n^{-1}\right)$, or alternatively $\hat{a}=\frac{n}{n+r+1}$. I present some numerical results below.
Placing last is finishing after $n=r$ clones finish:
[ n=1\longrightarrow\mathcal{L}=b]
[ n=2\longrightarrow\mathcal{L}=\frac{b}{a+2b}\cdot\frac{b}{a+b}=\frac{b^2}{1+b} =\frac{1}{B(B+1)} ]
[ n=3\longrightarrow\mathcal{L}=\frac{b}{a+3b}\cdot\frac{b}{a+2b}\cdot\frac{b}{a+b}=\frac{b^3}{(1+2b)(1+b)}=\frac{1}{B(B+1)(B+2)} ]
[ n=4\longrightarrow\mathcal{L}=\frac{b}{a+4b}\cdot\frac{b}{a+3b}\cdot\frac{b}{a+2b}\cdot\frac{b}{a+b}= \frac{b^4}{(1+3b)(1+2b)(1+b)}= \frac{1}{B(B+1)(B+2)(B+3)} ]
For $n=4$ we have equivalently
[ \frac{1}{(B+3)(B+2)(B+1)B}=\frac{(B-1)!}{(B+3)!} ]
In general we would have
[ \frac{1}{(B+n-1)(B+n-2)\cdots(B+1)B}= \frac{(B-1)!}{(B+n-1)!} ]
(denominator has $n+1$ terms). Observe that this agrees with the $r<n$ case above (substituting $r=n$). We may introduce some R idiom:
f_single <- function(a,r,n,log=FALSE){ # not vectorised B <- 1/(1-a) if(log){ out <- log(B-1) - lgamma(B+n) + lgamma(B+n-r-1) } else { # not-log out <- (B-1)/prod(B+(n-r-1):(n-1)) } return(out) } f_vec_a <- function(a,r,n, log=FALSE){ # vectorised in 'a' but not in 'r' sapply(a,function(a){f_single(a,r=r,n=n,log=log)}) } f_vec_arn <- function(a,r,n,log=FALSE){ ## vectorized in 'a'; treats 'r' and 'n' as ## vectors of independent observations M <- cbind(r,n) if(log){ out <- 0 for(i in seq_len(nrow(M))){out <- out + f_vec_a(a,r=M[i,1],n=M[i,2],log=TRUE)} } else { out <- 1 for(i in seq_len(nrow(M))){out <- out * f_vec_a(a,r=M[i,1],n=M[i,2],log=FALSE)} } return(out) } fapprox <- function(a,r,n,log=FALSE){ if(log){ return(log(B-1)-B*(r+1)/n) } else { return(B-1 - (r+1)*exp(B)/n) } }
fdash <- function(a,r,n,log=TRUE){ B <- 1/(1-a) out <- (1/(B-1) - psigamma(B+n) + psigamma(B+n-r-1))*B^2 if(!log){out <- out * f_vec_arn(a,r,n,log=FALSE)} return(out) }
Now verify fdash()
:
a <- 0.34 d <- 1e-3 r <- 4 n <- 10 c( numerical = (f_single(a+d/2,r,n,log=TRUE)-f_single(a-d/2,r,n,log=TRUE))/d, analytical = fdash(a,r,n) ) c( numerical = (f_single(a+d/2,r,n,log=FALSE)-f_single(a-d/2,r,n,log=FALSE))/d, analytical = fdash(a,r,n,log=FALSE) )
Now use the derivative to find the maximum likelihood point
MLE <- function(r,n,give=FALSE){ d <- 1e-6 out <- optimize(f_vec_arn,c(d,1-d),r=r,n=n,log=TRUE,maximum=TRUE) if(!give){out <- out$maximum} return(out) } MLEa <- function(r,n,give=FALSE,small=1e-4){ out <- uniroot(f=function(a){fdash(a,r,n)},interval=c(small,1-small)) if(!give){out <- out$root} return(out) } showdiff <- function(x,y){c(way1=x,way2=y,diff=x-y,logratio=log(x/y))} rbind( showdiff(MLE(1,9),MLEa(1,9)), showdiff(MLE(3,9),MLEa(3,9)), showdiff(MLE(3,7),MLEa(3,7)) )
rbind( showdiff(MLE(10, 100), 100 / 111), showdiff(MLE(10, 500), 500 / 511), showdiff(MLE(10,1000),1000 /1011), showdiff(MLE(10,5000),5000 /5011) )
rbind( showdiff(MLEa(10, 100), 100 / 111), showdiff(MLEa(10, 500), 500 / 511), showdiff(MLEa(10,1000),1000 /1011), showdiff(MLEa(10,5000),5000 /5011) )
Now verification
a <- 0.23423434 b <- 1-a B <- 1/b n <- 9 r <- 3 c( way1 = b^3*a/((a+9*b) * (a+8*b) * (a+7*b) * (a+6*b)), way2 = (B-1)/prod(B+(5:8)), way3 = f_vec_arn(a,r,n,log=FALSE) )
out <- list() for(n in 1:10){ out[[n]] <- rep(NA,n+1) for(r in 0:n){ if(r==0){ jj <- 1 } else if(r==n){n jj <- 0 } else { jj <- MLE(r,n) } out[[n]][r+1] <- jj } } out
plotterp <- function(...){ plot(NA,xlim=c(1,10),ylim=c(0,1),type="n",xlab="n",ylab=expression(hat(a))) for(n in 1:10){ for(r in 0:n){ points(n,out[[n]][r+1],pch=16) } for(r in 0:n){ text(n+0.2,out[[n]][r+1],r,cex=0.5,col='gray') } } } plotterp() pdf(file="dotprobs.pdf") plotterp() dev.off()
The likelihood calculations above can be confusing becuase likelihood is not the same as probability.
Consider, for example, $n=2$ so we have 3 competitors.
probability of focal competitor coming first [i.e. $r=0$]:
[\frac{a}{a+2b}=\frac{1-b}{1+b}]
Now what is the probability of coming second, that is, $r=1$?
My first thought was:
[ \mbox{prob of coming second} = \frac{b}{a+2b}\cdot\frac{a}{a+b}=\frac{b(1-b)}{1+b}\qquad\mbox{WRONG} ]
but this is incorrect: there are two clones of (strength) $b$, and the likelihood arguments above do not distinguish between their finishing order. We can name the competitors $a$, $b_1$ and $b_2$ and specify that $b_1$ and $b_2$ both have strength $b$. Then "coming second" means that the finishing order was $b_1\succ a\succ b_2$ or $b_2\succ a\succ b_1$. Noting that these events are disjoint, the probability would be
[\mbox{prob of coming second} = \operatorname{P}(b_1\succ a\succ b_2) + \operatorname{P}(b_2\succ a\succ b_1) = \frac{b}{a+2b}\cdot\frac{a}{a+b} + \frac{b}{a+2b}\cdot\frac{a}{a+b} = \frac{2(1-b)}{1+b} ]
Similarly for the focal competitor coming last we need to sum the probabilities of finishing orders $b_1\succ b_2\succ a$ and $b_2\succ b_1\succ a$:
[\mbox{prob of coming third}= \operatorname{P}(b_1\succ b_2\succ a) + \operatorname{P}(b_2\succ b_1\succ a) = \frac{b}{a+2b}\cdot\frac{b}{a+b}+ \frac{b}{a+2b}\cdot\frac{b}{a+b}= \frac{2b^2}{1+b}. ]
Above we see a factor of two difference between the probabilities and the likelihoods presented previously. Note that there is no such factor for the probability of the focal competitor finishing first:
[ \mbox{prob of coming first} = \operatorname{P}(a\succ b_1\succ b_2) + \operatorname{P}(a\succ b_2\succ b_1) = \frac{b}{a+2b}\cdot\frac{b}{b+b} + \frac{b}{a+2b}\cdot\frac{b}{b+b} = \frac{1-b}{1+b} ]
Now we see that the probabilities sum to
[ \frac{1-b}{1+b} + \frac{2b^2}{1+b} + \frac{2(1-b)}{1+b}=1 ]
Now try $r=3$ [i.e. 4 competitors]. Taking things one step at a time:
[ \mbox{prob of coming first} = 3!\cdot\frac{a}{a+3b}\cdot\frac{b}{3b}\cdot\frac{b}{2b}\cdot\frac{b}{b}= \frac{1-b}{1+2b} ]
[ \mbox{prob of coming second} = 3!\cdot\frac{b}{a+3b}\cdot\frac{a}{a+2b}\cdot\frac{b}{2b}\cdot\frac{b}{b}= 3\cdot\frac{b(1-b)}{(1+2b)(1+b)} ]
[ \mbox{prob of coming third} = 3!\cdot\frac{b}{a+3b}\cdot\frac{b}{a+2b}\cdot\frac{a}{a+b}\cdot\frac{b}{b} = 6\cdot\frac{b^2(1-b)}{(1+2b)(1+b)} ] [ \mbox{prob of coming fourth (last)} = 3!\cdot\frac{b}{a+3b}\cdot\frac{b}{a+2b}\cdot\frac{b}{a+b}\cdot\frac{a}{a} = 6\cdot\frac{b^3}{(1+2b)(1+b)} ]
The sum would be $[(1+b)(1-b) + 3b(1-b) + 6b^2(1-b) + 6b^3]/[(1+b)(1+2b)]=1$.
Now consider the general case. Again, $n$ clones of strength $b$ and one of strength $a=1-b$, $n+1$ competitors.
[ \mbox{prob of coming first} = n!\cdot\frac{b}{a+nb}\cdot\frac{b}{nb}\cdot\frac{b}{(n-1)b}\cdots\frac{b}{b} = \frac{b}{1+(n-1)b} ]
[ \mbox{prob of coming second} = n!\cdot\frac{b}{a+nb}\cdot\frac{a}{a+(n-1)b}\cdot\frac{b}{(n-1)b}\cdot\frac{b}{(n-2)b}\cdots\frac{b}{2b}\cdot\frac{b}{b} =n\cdot X ]
[ \mbox{prob of coming third} = \frac{b}{a+nb}\cdot \frac{b}{a+(n-1)b}\cdot \frac{a}{a+(n-2)b}\cdot \frac{b}{(n-2)b}\cdots \frac{b}{2b}\cdot\frac{b}{b} =\frac{n!}{(n-2)!}\cdot X ]
Now consider the case where there are $r$ clones coming ahead of the focal competitor:
[ \mbox{prob of having $r$ clones ahead} =\frac{n!}{(n-r)!}\cdot X ]
Now some likelihood visuals. Suppose $r=4,n=8$:
plotter <- function(...){ a <- seq(from=0,to=1,by=0.01) n <- 8 jj <- f_vec_arn(a,3,n,log=FALSE) plot(a,jj/max(jj,na.rm=TRUE),type='n',xlab=expression(a),ylab="likelihood") grid() rain <- rainbow(n+1) for(r in seq(from=0,to=n)){ y <- f_vec_a(a,r,n,log=FALSE) if(r==0){y[length(y)] <- 1} y <- y/max(y,na.rm=TRUE) if(r>0){y[length(y)] <- 0} points(a,y,type='l',lwd=4,col=rain[r+1]) } abline(v=MLE(r=4,n=8)) text(0.07,0.95,"r=8",col=rain[9]) text(0.20,0.95,"r=7",col=rain[8]) text(0.26,0.91,"r=6",col=rain[7]) text(0.29,0.83,"r=5",col=rain[6]) text(0.33,0.76,"r=4",col=rain[5]) text(0.35,0.64,"r=3",col=rain[4]) text(0.39,0.53,"r=2",col=rain[3]) text(0.45,0.41,"r=1",col=rain[2]) text(0.63,0.22,"r=0",col=rain[1]) } plotter() pdf(file="ninelikes.pdf") plotter() dev.off()
MLE(r=4,n=8)
n <- 7 a <- seq(from=0,to=1,by=0.01) fhyper3 <- function(r,n){ out <- hyper3() out['a'] <- 1 out['b'] <- r for(i in (n-r):n){ out[c(a=1,b=i)] %<>% dec } return(out) } plot(NA,xlim=c(0,1),ylim=c(0,1),type='n') M <- cbind(a=a,b=1-a) for(r in 0:7){ y <- loglik(M,fhyper3(r,n),log=FALSE) y <- y/max(y, na.rm=TRUE) points(M[,1],y,type="l",lwd=8) }
a <- read.table("formula1_2023.txt",header=TRUE) a <- a[,seq_len(ncol(a)-1)] d_ocon <- as.numeric(as.matrix(a)["Ocon",]) d_gasly <- as.numeric(as.matrix(a)["Gasly",]) d_ocon [is.na(d_ocon )] <- 22 d_gasly[is.na(d_gasly)] <- 22 d_ocon d_gasly MLE(d_ocon ,22) MLE(d_gasly,22)
two sample test:
both <- c(d_ocon,d_gasly) bothmax <- MLE(both,22) d_ocon_max <- MLE(d_ocon ,22) d_gasly_max <- MLE(d_gasly,22) f_vec_arn(bothmax,d_ocon ,22,log=TRUE) f_vec_arn(bothmax,d_gasly,22,log=TRUE) f_vec_arn(d_ocon_max ,d_ocon ,22,log=TRUE) f_vec_arn(d_gasly_max,d_gasly,22,log=TRUE) Lambda <- ( (f_vec_arn(d_gasly_max ,d_ocon,22,log=TRUE)+f_vec_arn(d_gasly_max ,d_gasly,22,log=TRUE)) - (f_vec_arn(bothmax ,d_ocon,22,log=TRUE)+f_vec_arn(bothmax ,d_gasly,22,log=TRUE)) ) Lambda pchisq(-2*Lambda,df=1,lower.tail=FALSE)
readstring <- function(year){read.table(paste("formula1_",year,".txt",sep=""))} getfoc <- function(year,comp="Perez"){ # get focal competitor M <- as.matrix(readstring(year)) out <- suppressWarnings(as.numeric(M[comp,seq_len(ncol(M)-1)])) out[is.na(out)] <- nrow(M) return(out) } getnum <- function(year){nrow(readstring(year))} # number of competitors perez <- lapply(2011:2023,getfoc) print(perez) y <- unlist(lapply(seq_along(perez),function(i){mean(perez[[i]])})) x <- 2011:2023 plot(x,y) summary(lm(y~x))
checo_like <- function(a){ out <- a*0 for(year in 2011:2023){ out <- out + f_vec_arn(a,getfoc(year)-1,getnum(year)-1,log=TRUE) } return(out) } a <- seq(from=0.45,to=0.62,by=0.01) cL <- checo_like(a) cL <- cL - max(cL) plot(a,cL,type='b') abline(h=c(0,-2))
f1_logistic <- function(vec){ alpha <- vec[1] beta <- vec[2] out <- 0 for(year in 2011:2023){ LO <- alpha + beta*(year-2011) strength <- exp(LO)/(1+exp(LO)) out <- out + f_vec_arn(strength,getfoc(year)-1,getnum(year)-1,log=TRUE) } return(out) }
a <- seq(from=-0.8,to=0.0,by=0.01) b <- seq(from=0,to=0.2,by=0.01) jj <- as.matrix(expand.grid(a,b)) L <- apply(jj,1,f1_logistic) L <- L-max(L) L <- matrix(L,length(a),length(b))
L <- pmax(L,-40) showchec <- function(...){ contour(a,b,L,xlab=expression(alpha),ylab=expression(beta),levels=-c(2,5*(1:5))) points(-0.27,0.0813,pch=16,cex=2) } showchec() pdf(file="showchecolike.pdf") showchec(a,b,L) dev.off()
f1_logistic(c(0,0)) f1_logistic(c(0,0.0001)) f1_logistic(c(0.001,0))
o <- optim(c(-0.2,0.1),fn=f1_logistic, control=list(fnscale = -1),hessian=TRUE)
o o$par o$hessian eigen(o$hessian) best <- o$par f1_logistic(best) jj <- best jj[1] <- jj[1]*1.01 f1_logistic(jj) - f1_logistic(best) jj <- best jj[1] <- jj[1]*0.99 f1_logistic(jj) - f1_logistic(best) jj <- best jj[2] <- jj[2]*1.00001 f1_logistic(jj) - f1_logistic(best) jj <- best jj[2] <- jj[2]*0.99 f1_logistic(jj) - f1_logistic(best)
o_free <- optim(c(-0.2,0.1),fn=f1_logistic, control=list(fnscale = -1),hessian=TRUE) o_constrained <- optim(c(-0.2,0.1),fn=function(v){ v[2] <- 0 return(f1_logistic(v))}, control=list(fnscale = -1),hessian=TRUE)
o_free o_constrained pchisq(2*(o_free$value - o_constrained$value),df=1,lower.tail=FALSE)
table(replicate(100,which(rrace3(pn=c(a=1,b=10),ps=c(a=0.9,b=0.1))=='a')))
n <- 1000 M <- matrix(0,n,4) alpha_true <- 4 beta_true <- -6 for(i in seq_len(n)){ x <- i/n LO <- alpha_true + beta_true*x p <- exp(LO)/(1+exp(LO)) nn <- round(runif(1,10,100)) M[i,1] <- x # regressor M[i,2] <- which(rrace3(pn=c(a=1,b=nn),ps=c(a=p,b=1-p))=='a') M[i,3] <- nn+1 M[i,4] <- p } head(M) tail(M)
likeinsilico <- function(vec,M){ alpha <- vec[1] beta <- vec[2] S <- 0 for(i in seq_len(nrow(M))){ x <- M[i,1] LO <- alpha + beta * x p <- exp(LO)/(1+exp(LO)) S <- S + f_vec_arn(p,M[i,2]-1,M[i,3]-1,log=TRUE) } return(S) }
likeinsilico(c(4,-6),M) likeinsilico(c(4,-6.5),M) likeinsilico(c(4,-5.5),M)
optim_ab <- optim(par=c(4,-6),likeinsilico,control=list(fnscale= -1),M=M) optim_ab
a <- seq(from=3.8,to=4.2,by=0.01) b <- seq(from=-6.4,to=-5.6,by=0.01) Z <- as.matrix(expand.grid(a,b)) Z <- apply(Z,1,function(x){likeinsilico(x,M=M)}) Z <- matrix(Z,length(a),length(b))
Z <- Z-max(Z) Z <- pmax(Z,-10) contour(a,b,Z,levels=-c(0:10,15,20,30)) points(alpha_true,beta_true,pch=16,cex=4) jj <- optim_ab$par points(jj[1],jj[2],pch=4,cex=4,col='red') grid() persp(a,b,Z)
sametest <- function(d1,d2,n){ both <- c(d1,d2) bothmax <- MLE(both,n) d1max <- MLE(d1,n) d2max <- MLE(d2,n) Lambda <- ( (f_vec_arn(d2max ,d1,n,log=TRUE)+f_vec_arn(d2max ,d2,n,log=TRUE)) - (f_vec_arn(bothmax,d1,n,log=TRUE)+f_vec_arn(bothmax,d2,n,log=TRUE)) ) pchisq(-2*Lambda,df=1,lower.tail=FALSE) } a <- as.matrix(read.table("olympiad.txt")) nrow(a) d1 <- a["AUS",] d2 <- a["NZL",] d1 d2 sametest(d1,d2,100) a["USA",] a["CHN",] sametest(a["USA",],a["CHN",],100) sametest(a["USA",],a["CHN",],5)
O <- read.table("olympic_athletics_mens_100m.txt",header=TRUE) head(O)
jjf <- function(x){f_vec_arn(x,O$rank,O$n,log=TRUE)} system.time(optolym <- optimize(jjf, c(0.45,0.5), maximum=TRUE))
optolym jjf(0.5) jjf(optolym$maximum) (LR <- 2*(jjf(optolym$maximum)-jjf(0.5))) pchisq(LR,df=1,lower.tail=FALSE)
a <- seq(from=0.3,to=0.65,len=45) L <- f_vec_arn(a,O$rank,O$n,log=TRUE) Lmax <- jjf(optolym$maximum) Lmax L <- L-Lmax
(a_lower <- uniroot(function(x){jjf(x)+2-Lmax},interval=c(0.3,0.5))$root) (a_upper <- uniroot(function(x){jjf(x)+2-Lmax},interval=c(0.5,0.6))$root)
plotolymp <- function(...){ plot(a,L,type='l',ylab="log-likelihood",lwd=2) abline(h=c(0,-2)) ahat <- optolym$maximum segments(x0=ahat,y0=-1.5,y1=0) text(ahat,-0.91,expression(hat(a)==0.573),pos=2) abline(v=0.5,lty=3) segments(x0=a_lower,y0=-1.5,y1=-3.5) segments(x0=a_upper,y0=-1.5,y1=-3.5) text(x=a_lower,y=-3,"0.380",pos=4) text(x=a_upper,y=-3,"0.563",pos=2) } plotolymp() pdf(file="plotolymp.pdf") plotolymp() dev.off()
results <- read.table("parkrun_results.txt",header=TRUE) results
Above we see some parkrun results. In parkrun number 148, I placed 229 out of a field of 318.
S_parkrun <- function(a){f_vec_arn(a,r=results$place-1,results$runners-1,log=TRUE)} (optparkrun <- optimize(S_parkrun, c(0.2,0.8),maximum=TRUE))
a <- seq(from=0.25,to=0.65,by=0.01) S_parkrun_a <- S_parkrun(a) - optparkrun$objective uniroot(function(a){S_parkrun(a) - optparkrun$objective+2},c(0.1,0.5))$root uniroot(function(a){S_parkrun(a) - optparkrun$objective+2},c(0.5,0.9))$root optparkrun$objective - S_parkrun(0.5)
```r\simeq 0.448$ has a support of zero. Dotted line shows $\mathcal{S}=-2$; intersection of this with the curve shows that the support interval runs from about 0.315 to 0.567"} plotparkrun <- function(...){ plot(a,S_parkrun_a,xlab="generalized Bradley-Terry strength",ylab="support",type='b',ylim=c(-5,0)) abline(h=c(0,-2),lty=2) segments(x0=0.5,y0=-1.5,y1=0) segments(x0=optparkrun$maximum,y0=-1,y1=0) text(optparkrun$maximum,-0.8,expression(hat(a)==0.448),pos=2) text(0.5,-1.5,expression(H[o]:a==0.5),pos=2) } plotparkrun() pdf("plotparkrun.pdf") plotparkrun() dev.off()
Figure \@ref(fig:plotparkrunlike) shows a support curve for the author's generalized Bradley-Terry strength, normalized so that $\mathcal{S}(\hat{a})=0$. We see a support interval [that is, $\left\lbrace a\colon\mathcal{S}(a)\geqslant -2\right\rbrace$] of about 0.315 to 0.567. We might also observe that the support for $H_0\colon a=\frac{1}{2}$ [that is, that the author has a 50\% chance of running faster than a randomly chosen Parkrun competitor] is about 0.35, less than Edwards's two units of support criterion. We thus see some support for Hankin's assertion that he is "as fit as anyone else there". ```r a <- seq(from=0.005,to=0.85,by=0.01) La <- f_vec_arn(a,r=results$place[1]-1,results$runners[1]-1,log=TRUE) La <- La - max(La,na.rm=TRUE) plot(a,La,type='b',ylab='fishy') abline(h=c(0,-2))
supp2 <- function(vec,place,runners){ a_stirling <- vec[1] a_auckland <- vec[2] out <- 0 for(i in 1:3){ out <- out + f_vec_a(a_stirling,r=place[i]-1,runners[i]-1,log=TRUE) } for(i in 4:10){ out <- out + f_vec_a(a_auckland,r=place[i]-1,runners[i]-1,log=TRUE) } return(out) }
x <- seq(from=0.1,to=0.8,by=0.01) jj <- as.matrix(expand.grid(x,x)) L <- apply(jj,1,supp2,results$place,results$runners) L <- L - max(L,na.rm=TRUE) L <- matrix(L,length(x),length(x))
par(pty='s') contour(x,x,L,asp=1,xlim=range(x),ylim=range(x),levels=-(0:5)) abline(0,1)
f <- function(a){supp2(c(a,a),results$place,results$runners)} a <- seq(from=0.2,to=0.6,by=0.01) alike <- sapply(a,f) plot(a,alike-max(alike,na.rm=TRUE),type='b',ylab='blong') abline(h=c(0,-2)) optimize(f,c(0.1,0.6),maximum=TRUE)
f2 <- function(v){supp2(v,results$place,results$runners)} optim(par=c(0.6,0.3),f2,control=list(fnscale=-1))
7380.314-7381.496 pchisq(2*(7381.496-7380.314),df=1,lower.tail=FALSE)
Above shows that there is no evidence to suggest that the Stirling strength differs from the Auckland strength. Now let's try a logistic-style regression:
results <- read.table("parkrun_results.txt",header=TRUE) results supp_logistic <- function(vec){ alpha <- vec[1] beta <- vec[2] out <- 0 for(i in seq_len(nrow(results))){ LO <- alpha + beta*i strength <- exp(LO)/(1+exp(LO)) out <- out + f_vec_arn(strength,results$place[i]-1,results$runners[i]-1,log=TRUE) } return(out+13400) } c( supp_logistic(c(0,0)), supp_logistic(c(0,0.01)), supp_logistic(c(0,-0.01)) )
optim(par=c(0,0),fn=supp_logistic,control=list(fnscale=-1)) optim(par=c(0,0),fn=function(v){supp_logistic(vec=c(v[1],0))},control=list(fnscale=-1))
The focal competitor
rank <- c( 9, 7, 2, 3,2 ,8 ,5 ,4 ,9 ) class_size <- c(12, 17,23, 9,13,14,13,12,15) course <- c("rings and modules","group theory","calculus","linear algebra", "differential equations","topology","special relativity","fluid mechanics","Lie algebra") category=c("pure","pure","applied","applied","applied","pure","applied","applied","pure") data.frame(course,category,rank,class_size)
a <- seq(from=0.2,by=0.01,to=0.8) wp <- category=='pure' wa <- category=='applied'
Pure first
Lpure <- f_vec_arn(a,rank[wp]-1,class_size[wp]-1,log=TRUE) Lpure <- Lpure - max(Lpure,na.rm=TRUE) plot(a,Lpure,ylab='pure strength') optimize(function(a){f_vec_arn(a,rank[wp]-1,class_size[wp]-1,log=TRUE)},c(0.3,0.6),maximum=TRUE) optimize(function(a){f_vec_arn(a,rank[wa]-1,class_size[wa]-1,log=TRUE)},c(0.3,0.6),maximum=TRUE)
supp2_ed <- function(vec,place1,place2,runners1,runners2){ out <- 0 M1 <- cbind(place1,runners1) M2 <- cbind(place2,runners2) for(i in seq_len(nrow(M1))){ out <- out + f_vec_arn(vec[1],r=M1[i,1]-1,n=M1[i,2]-1,log=TRUE) } for(i in seq_len(nrow(M2))){ out <- out + f_vec_arn(vec[2],r=M2[i,1]-1,n=M2[i,2]-1,log=TRUE) } return(out) }
x <- seq(from=0.1,to=0.99,by=0.01) jj <- as.matrix(expand.grid(x,x)) L <- apply(jj,1,supp2_ed,rank[wp]-1,rank[wa]-1,class_size[wp]-1,class_size[wa]-1) L <- L - max(L,na.rm=TRUE) L <- matrix(L,length(x),length(x))
plotpureandapplied <- function(...){ par(pty='s') contour(x,x,L,asp=1,xlim=range(x),ylim=range(x),levels=-(0:5), xlab='pure strength',ylab='applied strength') abline(0,1) } plotpureandapplied() pdf(file="plotpureandapplied.pdf") plotpureandapplied() dev.off()
jj1 <- optim( c(.4,.8), function(vec){ supp2_ed(vec,rank[wp]-1,rank[wa]-1,class_size[wp]-1,class_size[wa]-1)}, control=list(fnscale = -1) ) jj1
jj2 <- optimize(function(v){supp2_ed(c(v,v),rank[wp]-1,rank[wa]-1,class_size[wp]-1,class_size[wa]-1)}, c(0.1,0.9),maximum=TRUE) jj2
pchisq(2*(jj1$value-jj2$objective),lower.tail=FALSE,df=1)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.