(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)
##        1        2        3        4 
## 0.533333 0.266667 0.133333 0.066667
sum(BT(4,0.5))
## [1] 1

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)
##       [,1] [,2] [,3] [,4] [,5]
##  [1,]    2    1    3    4    5
##  [2,]    2    3    4    5    1
##  [3,]    4    5    1    3    2
##  [4,]    2    3    1    4    5
##  [5,]    3    2    1    5    4
##  [6,]    1    2    4    3    5
##  [7,]    2    3    5    1    4
##  [8,]    1    3    2    4    5
##  [9,]    1    2    3    4    5
## [10,]    1    2    3    4    5
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)
##  [1]  996 1193 1347 1468 1566 1342 1065  620  293   89   21
get_one_competitor(X,11)
##  [1]   74  110  140  195  233  358  519  754 1230 2055 4332
rbind(
get_one_competitor(X,1) ,
get_one_competitor(X,2) ,
get_one_competitor(X,9) ,
get_one_competitor(X,10)
)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 3024 2457 1814 1323  753  404  169   45    9     2     0
## [2,] 2162 2142 1887 1484 1082  720  353  127   35     8     0
## [3,]  175  206  270  335  486  698  948 1429 1778  2129  1546
## [4,]  130  135  180  283  316  469  720 1106 1607  2377  2677
t(sapply(seq_len(11),function(i){get_one_competitor(X,i)}))
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
##  [1,] 3024 2457 1814 1323  753  404  169   45    9     2     0
##  [2,] 2162 2142 1887 1484 1082  720  353  127   35     8     0
##  [3,] 1524 1582 1666 1629 1421 1037  663  339  120    14     5
##  [4,]  996 1193 1347 1468 1566 1342 1065  620  293    89    21
##  [5,]  758  874 1035 1172 1423 1495 1334  993  606   272    38
##  [6,]  542  594  722  947 1183 1420 1529 1399  991   531   142
##  [7,]  367  425  540  682  866 1152 1453 1621 1549   937   408
##  [8,]  248  282  399  482  671  905 1247 1567 1782  1586   831
##  [9,]  175  206  270  335  486  698  948 1429 1778  2129  1546
## [10,]  130  135  180  283  316  469  720 1106 1607  2377  2677
## [11,]   74  110  140  195  233  358  519  754 1230  2055  4332

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)
##          observed_rank
## true_rank  c_1  c_2  c_3  c_4  c_5  c_6  c_7  c_8  c_9
##       r_1 0.54 0.28 0.10 0.06 0.02 0.00 0.00 0.00 0.00
##       r_2 0.22 0.32 0.27 0.15 0.02 0.02 0.00 0.00 0.00
##       r_3 0.16 0.22 0.33 0.15 0.08 0.06 0.00 0.00 0.00
##       r_4 0.05 0.07 0.11 0.29 0.30 0.14 0.04 0.00 0.00
##       r_5 0.01 0.05 0.09 0.20 0.25 0.23 0.15 0.02 0.00
##       r_6 0.01 0.04 0.05 0.09 0.15 0.27 0.21 0.15 0.03
##       r_7 0.01 0.02 0.03 0.03 0.10 0.12 0.34 0.24 0.11
##       r_8 0.00 0.00 0.02 0.02 0.05 0.06 0.15 0.40 0.30
##       r_9 0.00 0.00 0.00 0.01 0.03 0.10 0.11 0.19 0.56

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)
##          observed_rank
## true_rank  c_1  c_2  c_3  c_4  c_5  c_6  c_7  c_8  c_9
##       r_1 0.11 0.14 0.12 0.08 0.10 0.06 0.17 0.12 0.10
##       r_2 0.07 0.05 0.10 0.14 0.19 0.08 0.16 0.10 0.11
##       r_3 0.13 0.10 0.11 0.11 0.15 0.12 0.09 0.11 0.08
##       r_4 0.14 0.15 0.13 0.15 0.06 0.14 0.06 0.10 0.07
##       r_5 0.15 0.12 0.11 0.10 0.11 0.13 0.10 0.10 0.08
##       r_6 0.06 0.12 0.17 0.07 0.09 0.12 0.09 0.12 0.16
##       r_7 0.14 0.12 0.10 0.13 0.12 0.10 0.13 0.10 0.06
##       r_8 0.15 0.08 0.10 0.11 0.09 0.11 0.08 0.15 0.13
##       r_9 0.05 0.12 0.06 0.11 0.09 0.14 0.12 0.10 0.21

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

get_all_competitors(100, N, 0.1)
##          observed_rank
## true_rank  c_1  c_2  c_3  c_4  c_5  c_6  c_7  c_8 c_9
##       r_1 0.89 0.11 0.00 0.00 0.00 0.00 0.00 0.00 0.0
##       r_2 0.10 0.81 0.09 0.00 0.00 0.00 0.00 0.00 0.0
##       r_3 0.01 0.07 0.83 0.09 0.00 0.00 0.00 0.00 0.0
##       r_4 0.00 0.01 0.08 0.88 0.03 0.00 0.00 0.00 0.0
##       r_5 0.00 0.00 0.00 0.02 0.88 0.10 0.00 0.00 0.0
##       r_6 0.00 0.00 0.00 0.01 0.09 0.82 0.08 0.00 0.0
##       r_7 0.00 0.00 0.00 0.00 0.00 0.08 0.86 0.06 0.0
##       r_8 0.00 0.00 0.00 0.00 0.00 0.00 0.05 0.85 0.1
##       r_9 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.09 0.9

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)
##        x r     c_1     c_2     c_3     c_4     c_5     c_6     c_7     c_8 c_9
## [1,] 0.1 1 0.90074 0.09728 0.00198 0.00000 0.00000 0.00000 0.00000 0.00000   0
## [2,] 0.1 2 0.08878 0.81304 0.09618 0.00200 0.00000 0.00000 0.00000 0.00000   0
## [3,] 0.1 3 0.00950 0.08106 0.80994 0.09754 0.00196 0.00000 0.00000 0.00000   0
## [4,] 0.1 4 0.00092 0.00794 0.08242 0.80608 0.10068 0.00196 0.00000 0.00000   0
## [5,] 0.1 5 0.00006 0.00054 0.00880 0.08488 0.80810 0.09550 0.00210 0.00002   0
## [6,] 0.1 6 0.00000 0.00010 0.00058 0.00834 0.08042 0.81260 0.09618 0.00178   0
allX9[seq(from=round(nrow(allX9)/2 - 5), to=round(nrow(allX9)/2 + 5)),]
##         x r     c_1     c_2     c_3     c_4     c_5     c_6     c_7     c_8
##  [1,] 0.4 9 0.00042 0.00080 0.00182 0.00480 0.01128 0.02862 0.07210 0.20588
##  [2,] 0.5 1 0.50060 0.30346 0.14198 0.04484 0.00808 0.00100 0.00004 0.00000
##  [3,] 0.5 2 0.24864 0.31636 0.24878 0.13136 0.04468 0.00900 0.00110 0.00008
##  [4,] 0.5 3 0.12436 0.18400 0.27526 0.23980 0.12584 0.04154 0.00840 0.00072
##  [5,] 0.5 4 0.06446 0.09926 0.16210 0.26872 0.23346 0.12576 0.03926 0.00640
##  [6,] 0.5 5 0.03214 0.05024 0.08794 0.15636 0.27604 0.23510 0.12258 0.03546
##  [7,] 0.5 6 0.01554 0.02410 0.04330 0.08388 0.16212 0.28348 0.24502 0.11848
##  [8,] 0.5 7 0.00836 0.01326 0.02320 0.04338 0.08384 0.16804 0.30860 0.25346
##  [9,] 0.5 8 0.00384 0.00614 0.01128 0.02108 0.04444 0.08944 0.17762 0.37124
## [10,] 0.5 9 0.00206 0.00318 0.00616 0.01058 0.02150 0.04664 0.09738 0.21416
## [11,] 0.6 1 0.40768 0.28542 0.17568 0.08762 0.03344 0.00842 0.00166 0.00008
##           c_9
##  [1,] 0.67428
##  [2,] 0.00000
##  [3,] 0.00000
##  [4,] 0.00008
##  [5,] 0.00058
##  [6,] 0.00414
##  [7,] 0.02408
##  [8,] 0.09786
##  [9,] 0.27492
## [10,] 0.59834
## [11,] 0.00000
tail(allX9)
##         x r     c_1     c_2     c_3     c_4     c_5     c_6     c_7     c_8
## [76,] 0.9 4 0.11956 0.11910 0.11852 0.11784 0.11810 0.11430 0.10996 0.10190
## [77,] 0.9 5 0.10620 0.10882 0.11074 0.11336 0.11444 0.11870 0.11542 0.11344
## [78,] 0.9 6 0.09622 0.10132 0.10166 0.10630 0.11248 0.11332 0.12034 0.12666
## [79,] 0.9 7 0.08424 0.09046 0.09568 0.09990 0.10284 0.11450 0.12360 0.13574
## [80,] 0.9 8 0.07954 0.08044 0.08700 0.09232 0.10018 0.10944 0.12446 0.14532
## [81,] 0.9 9 0.07068 0.07476 0.07942 0.08736 0.09506 0.10380 0.12218 0.14972
##           c_9
## [76,] 0.08072
## [77,] 0.09888
## [78,] 0.12170
## [79,] 0.15304
## [80,] 0.18130
## [81,] 0.21702

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)
## [1] 2 0 3 2 1 4 2 3 3
make_obs(no_of_observations=20, no_of_competitors=N, x_true = 0.7, r_true = 9)
## [1]  0  1  0  0  0  1  1  4 13
make_obs(no_of_observations=20, no_of_competitors=N, x_true = 0.1, r_true = 6)
## [1]  0  0  0  0  2 16  2  0  0
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))
##      r
## x      1  2  3  4  5  6  7  8  9
##   0.3  0  0  1  0  0  0  0  0  0
##   0.4  0  0  2  4  6  4  0  0  0
##   0.5  0  0  5 18 29 12  1  0  0
##   0.6  1  3 13 59 67 14  5  0  0
##   0.7  3  9 47 73 68 39  2  1  0
##   0.8 23 36 64 46 39 23  5  1  0
##   0.9 54 38 48 45 25 37 19  6  5
biasx(jj, x_true = 0.8)
## [1] -0.0553
biasr(jj, r_true = 4)
## [1] 0.05
msex(jj, x_true = 0.8)
## [1] 0.02033
mser(jj, r_true = 4)
## [1] 2.544

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
##      r
## x      1  2  3  4  5  6
##   0.3  0  0  0  0  1  1
##   0.4  0  0  1  3  0  0
##   0.5  0  0  0  6  5  2
##   0.6  0  2  3 12  5  4
##   0.7  0  1  6 11  4  1
##   0.8  1  5  3  4  4  0
##   0.9  5  3  3  1  2  1
jj2
##      r
## x      1  2  3  4  5  6  7
##   0.4  0  0  0  1  2  2  0
##   0.5  0  0  2  6  2  1  1
##   0.6  0  0  5 11 12  2  0
##   0.7  0  2  1 10  3  2  1
##   0.8  4  5  8  5  3  1  0
##   0.9  3  1  1  0  3  0  0
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
##       bias(r) mse(r) bias(x) mse(x)
##  [1,]   -0.17  1.682  -0.025 0.0223
##  [2,]   -0.07  1.882  -0.027 0.0181
##  [3,]   -0.14  1.644  -0.024 0.0190
##  [4,]   -0.16  1.776  -0.026 0.0196
##  [5,]   -0.12  1.472  -0.039 0.0177
##  [6,]   -0.05  1.870  -0.012 0.0200
##  [7,]   -0.07  1.502  -0.052 0.0222
##  [8,]   -0.34  2.004  -0.019 0.0171
##  [9,]   -0.24  1.624  -0.018 0.0198
mean(fish[,1])
## [1] -0.15111
mean(fish[,2])
## [1] 1.7173
sd(fish[,1])
## [1] 0.092526
mean(fish[,2])/sqrt(100)
## [1] 0.17173
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
##      r
## x       1   2   3   4   5   6   7
##   0.2   0   0   0   1   0   0   0
##   0.3   0   0   2   0   2   0   0
##   0.4   0   1   1  19  12   5   0
##   0.5   0   1  14  63  66  12   1
##   0.6   0   8  39 108  91   8   1
##   0.7   4  17  59 100  60  14   0
##   0.8  31  50  56  33  12   4   1
##   0.9  32  27  15  19   5   4   2
jj2
##      r
## x       1   2   3   4   5   6   7
##   0.2   0   0   0   0   0   1   0
##   0.3   0   0   0   3   6   2   0
##   0.4   0   0   6  13  15   5   1
##   0.5   0   2  15  56  49   9   0
##   0.6   0   8  37 103  91  13   0
##   0.7   3  17  63  90  67  18   0
##   0.8  27  35  57  32  17   9   0
##   0.9  44  25  23  18  12   7   1
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
##       bias(r) mse(r) bias(x)  mse(x)
##  [1,]  -0.238 1.6308 -0.0309 0.01917
##  [2,]  -0.206 1.7036 -0.0246 0.02034
##  [3,]  -0.194 1.6684 -0.0355 0.01983
##  [4,]  -0.213 1.6298 -0.0325 0.01995
##  [5,]  -0.217 1.6642 -0.0284 0.01850
##  [6,]  -0.211 1.7126 -0.0226 0.01824
##  [7,]  -0.296 1.7336 -0.0250 0.01930
##  [8,]  -0.157 1.6362 -0.0301 0.01937
##  [9,]  -0.147 1.6742 -0.0367 0.02087
mean(fish[,1])
## [1] -0.20878
mean(fish[,2])
## [1] 1.6726
sd(fish[,1])
## [1] 0.043697
mean(fish[,2])/sqrt(1000)
## [1] 0.052892

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
##      bias(r)   bias(x)  mse(r)    mse(x)
## 0.1 -0.00150  0.018500 0.00200 0.0022250
## 0.2 -0.01900 -0.004975 0.02750 0.0075275
## 0.3 -0.03875 -0.004250 0.10075 0.0100950
## 0.4 -0.04400 -0.013175 0.22150 0.0120625
## 0.5 -0.08975 -0.016275 0.47675 0.0148425
## 0.6 -0.16100 -0.022125 0.96700 0.0172475
## 0.7 -0.18325 -0.027325 1.76125 0.0197025
## 0.8  0.03550 -0.055675 2.84100 0.0204225
## 0.9  0.59125 -0.111525 4.47925 0.0276825
matplot(jj,type="b")

cbind(`bias(r)`=jj[,1], `SE(bias)` = jj[,3]/sqrt(4000))
##      bias(r)   SE(bias)
## 0.1 -0.00150 3.1623e-05
## 0.2 -0.01900 4.3481e-04
## 0.3 -0.03875 1.5930e-03
## 0.4 -0.04400 3.5022e-03
## 0.5 -0.08975 7.5381e-03
## 0.6 -0.16100 1.5290e-02
## 0.7 -0.18325 2.7848e-02
## 0.8  0.03550 4.4920e-02
## 0.9  0.59125 7.0823e-02
cbind(`bias(x)`=jj[,2], `SE(bias)` = jj[,4]/sqrt(4000))
##       bias(x)   SE(bias)
## 0.1  0.018500 0.00003518
## 0.2 -0.004975 0.00011902
## 0.3 -0.004250 0.00015962
## 0.4 -0.013175 0.00019072
## 0.5 -0.016275 0.00023468
## 0.6 -0.022125 0.00027271
## 0.7 -0.027325 0.00031152
## 0.8 -0.055675 0.00032291
## 0.9 -0.111525 0.00043770

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
##     bias(r) bias(x) mse(r)  mse(x)
## 0.1   0.634 -0.1058  0.970 0.02812
## 0.2   0.189 -0.0727  0.991 0.02171
## 0.3  -0.097 -0.0480  1.283 0.02094
## 0.4  -0.103 -0.0311  1.613 0.02125
## 0.5  -0.103 -0.0289  1.825 0.02053
## 0.6   0.070 -0.0258  1.744 0.01964
## 0.7   0.100 -0.0309  1.354 0.02003
## 0.8   0.011 -0.0466  0.809 0.02008
## 0.9  -0.367 -0.0721  0.535 0.02297
matplot(jj,type="b")

cbind(`bias(r)`=jj[,1], `SE(bias)` = jj[,3]/sqrt(1000))
##     bias(r) SE(bias)
## 0.1   0.634 0.030674
## 0.2   0.189 0.031338
## 0.3  -0.097 0.040572
## 0.4  -0.103 0.051008
## 0.5  -0.103 0.057712
## 0.6   0.070 0.055150
## 0.7   0.100 0.042817
## 0.8   0.011 0.025583
## 0.9  -0.367 0.016918
cbind(`bias(x)`=jj[,2], `SE(bias)` = jj[,4]/sqrt(1000))
##     bias(x)   SE(bias)
## 0.1 -0.1058 0.00088923
## 0.2 -0.0727 0.00068653
## 0.3 -0.0480 0.00066218
## 0.4 -0.0311 0.00067198
## 0.5 -0.0289 0.00064922
## 0.6 -0.0258 0.00062107
## 0.7 -0.0309 0.00063340
## 0.8 -0.0466 0.00063499
## 0.9 -0.0721 0.00072638

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
##     bias(r) bias(x) mse(r)  mse(x)
## 10   -0.196 -0.0280  1.714 0.01966
## 20   -0.286 -0.0020  1.260 0.01126
## 50   -0.192  0.0073  0.648 0.00585
## 100  -0.114  0.0058  0.308 0.00354
## 200  -0.071  0.0059  0.091 0.00143
matplot(jj,type="b")

cbind(`bias(r)`=jj[,1], `SE(bias)` = jj[,3]/sqrt(1000))
##     bias(r)  SE(bias)
## 10   -0.196 0.0542014
## 20   -0.286 0.0398447
## 50   -0.192 0.0204916
## 100  -0.114 0.0097398
## 200  -0.071 0.0028777
cbind(`bias(x)`=jj[,2], `SE(bias)` = jj[,4]/sqrt(1000))
##     bias(x)   SE(bias)
## 10  -0.0280 6.2170e-04
## 20  -0.0020 3.5607e-04
## 50   0.0073 1.8499e-04
## 100  0.0058 1.1194e-04
## 200  0.0059 4.5221e-05