(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.53333333 0.26666667 0.13333333 0.06666667
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,]    1    3    4    2    5
##  [2,]    3    2    1    4    5
##  [3,]    2    4    3    1    5
##  [4,]    2    1    3    5    4
##  [5,]    4    2    1    5    3
##  [6,]    5    2    3    1    4
##  [7,]    5    3    2    1    4
##  [8,]    5    2    1    3    4
##  [9,]    4    3    1    2    5
## [10,]    1    5    2    3    4
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] 1043 1234 1317 1477 1539 1325 1052  607  292   93   21
get_one_competitor(X,11)
##  [1]   84  116  142  189  236  332  511  774 1270 2024 4322
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,] 3084 2455 1867 1243  783  367  135   51   12     3     0
## [2,] 2110 2044 1899 1576 1117  706  378  127   38     4     1
## [3,]  192  214  254  364  506  688  976 1343 1806  2137  1520
## [4,]  125  145  189  261  357  511  717 1032 1579  2389  2695
t(sapply(seq_len(11),function(i){get_one_competitor(X,i)}))
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
##  [1,] 3084 2455 1867 1243  783  367  135   51   12     3     0
##  [2,] 2110 2044 1899 1576 1117  706  378  127   38     4     1
##  [3,] 1513 1595 1658 1619 1347 1122  664  345  107    28     2
##  [4,] 1043 1234 1317 1477 1539 1325 1052  607  292    93    21
##  [5,]  750  878 1032 1167 1427 1435 1310 1081  621   243    56
##  [6,]  537  578  760  969 1158 1410 1536 1354 1029   521   148
##  [7,]  336  435  508  639  883 1172 1487 1678 1464   998   400
##  [8,]  226  306  374  496  647  932 1234 1608 1782  1560   835
##  [9,]  192  214  254  364  506  688  976 1343 1806  2137  1520
## [10,]  125  145  189  261  357  511  717 1032 1579  2389  2695
## [11,]   84  116  142  189  236  332  511  774 1270  2024  4322

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.45 0.40 0.08 0.05 0.02 0.00 0.00 0.00 0.00
##       r_2 0.21 0.27 0.26 0.21 0.05 0.00 0.00 0.00 0.00
##       r_3 0.13 0.21 0.31 0.17 0.14 0.04 0.00 0.00 0.00
##       r_4 0.15 0.07 0.19 0.19 0.27 0.10 0.03 0.00 0.00
##       r_5 0.01 0.02 0.07 0.18 0.26 0.26 0.15 0.05 0.00
##       r_6 0.04 0.01 0.03 0.11 0.10 0.31 0.26 0.12 0.02
##       r_7 0.01 0.02 0.05 0.06 0.08 0.12 0.24 0.34 0.08
##       r_8 0.00 0.00 0.01 0.02 0.04 0.11 0.19 0.25 0.38
##       r_9 0.00 0.00 0.00 0.01 0.04 0.06 0.13 0.24 0.52

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.08 0.12 0.06 0.12 0.12 0.11 0.12 0.19 0.08
##       r_2 0.14 0.09 0.10 0.10 0.15 0.09 0.09 0.10 0.14
##       r_3 0.07 0.21 0.12 0.12 0.10 0.12 0.12 0.07 0.07
##       r_4 0.11 0.10 0.06 0.11 0.15 0.09 0.13 0.11 0.14
##       r_5 0.11 0.11 0.20 0.11 0.08 0.11 0.06 0.11 0.11
##       r_6 0.09 0.19 0.06 0.09 0.12 0.16 0.08 0.09 0.12
##       r_7 0.13 0.08 0.14 0.11 0.10 0.08 0.16 0.08 0.12
##       r_8 0.12 0.05 0.16 0.16 0.09 0.07 0.16 0.08 0.11
##       r_9 0.15 0.05 0.10 0.08 0.09 0.17 0.08 0.17 0.11

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.94 0.06 0.00 0.00 0.00 0.00 0.00 0.00 0.00
##       r_2 0.06 0.87 0.07 0.00 0.00 0.00 0.00 0.00 0.00
##       r_3 0.00 0.06 0.82 0.12 0.00 0.00 0.00 0.00 0.00
##       r_4 0.00 0.01 0.10 0.73 0.16 0.00 0.00 0.00 0.00
##       r_5 0.00 0.00 0.01 0.12 0.81 0.06 0.00 0.00 0.00
##       r_6 0.00 0.00 0.00 0.03 0.03 0.86 0.08 0.00 0.00
##       r_7 0.00 0.00 0.00 0.00 0.00 0.07 0.84 0.09 0.00
##       r_8 0.00 0.00 0.00 0.00 0.00 0.01 0.06 0.81 0.12
##       r_9 0.00 0.00 0.00 0.00 0.00 0.00 0.02 0.10 0.88

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.90124 0.09702 0.00174 0.00000 0.00000 0.00000 0.00000 0.0000   0
## [2,] 0.1 2 0.08894 0.81162 0.09738 0.00204 0.00002 0.00000 0.00000 0.0000   0
## [3,] 0.1 3 0.00884 0.08178 0.81122 0.09636 0.00180 0.00000 0.00000 0.0000   0
## [4,] 0.1 4 0.00090 0.00880 0.08042 0.81124 0.09672 0.00190 0.00002 0.0000   0
## [5,] 0.1 5 0.00008 0.00074 0.00846 0.08096 0.80868 0.09936 0.00172 0.0000   0
## [6,] 0.1 6 0.00000 0.00004 0.00068 0.00854 0.08398 0.80718 0.09748 0.0021   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.00038 0.00096 0.00184 0.00420 0.01096 0.02996 0.07370 0.20418
##  [2,] 0.5 1 0.50146 0.30136 0.14172 0.04500 0.00940 0.00100 0.00006 0.00000
##  [3,] 0.5 2 0.25110 0.32230 0.24746 0.12770 0.04264 0.00792 0.00084 0.00002
##  [4,] 0.5 3 0.12606 0.18200 0.27820 0.23474 0.12816 0.04178 0.00810 0.00094
##  [5,] 0.5 4 0.06262 0.09716 0.16200 0.27280 0.23368 0.12492 0.04048 0.00590
##  [6,] 0.5 5 0.03130 0.05000 0.08546 0.15966 0.27132 0.24198 0.12192 0.03420
##  [7,] 0.5 6 0.01432 0.02592 0.04512 0.08120 0.16196 0.28088 0.24986 0.11580
##  [8,] 0.5 7 0.00760 0.01218 0.02366 0.04456 0.08794 0.16934 0.30282 0.25704
##  [9,] 0.5 8 0.00390 0.00640 0.01088 0.02310 0.04426 0.08842 0.17816 0.36810
## [10,] 0.5 9 0.00164 0.00268 0.00550 0.01124 0.02064 0.04376 0.09776 0.21800
## [11,] 0.6 1 0.40348 0.28718 0.17580 0.08874 0.03326 0.00948 0.00178 0.00026
##           c_9
##  [1,] 0.67382
##  [2,] 0.00000
##  [3,] 0.00002
##  [4,] 0.00002
##  [5,] 0.00044
##  [6,] 0.00416
##  [7,] 0.02494
##  [8,] 0.09486
##  [9,] 0.27678
## [10,] 0.59878
## [11,] 0.00002
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.11880 0.11628 0.11840 0.11762 0.11604 0.11694 0.11068 0.10242
## [77,] 0.9 5 0.10830 0.11112 0.11136 0.11134 0.11326 0.11574 0.11762 0.11264
## [78,] 0.9 6 0.09874 0.09952 0.10198 0.10564 0.10916 0.11434 0.12068 0.12652
## [79,] 0.9 7 0.08416 0.08824 0.09544 0.10232 0.10772 0.11444 0.12054 0.13678
## [80,] 0.9 8 0.07698 0.08256 0.08914 0.09202 0.10012 0.10988 0.12384 0.14372
## [81,] 0.9 9 0.06984 0.07624 0.08028 0.08642 0.09406 0.10496 0.12166 0.15066
##           c_9
## [76,] 0.08282
## [77,] 0.09862
## [78,] 0.12342
## [79,] 0.15036
## [80,] 0.18174
## [81,] 0.21588

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] 1 2 3 1 3 4 2 3 1
make_obs(no_of_observations=20, no_of_competitors=N, x_true = 0.7, r_true = 9)
## [1]  0  0  0  0  1  2  1  4 12
make_obs(no_of_observations=20, no_of_competitors=N, x_true = 0.1, r_true = 6)
## [1]  0  0  0  1  3 15  1  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.4  0  0  2  6  8  5  0  0  0
##   0.5  0  0  3 20 19  9  1  0  0
##   0.6  1  3 22 53 66 23  5  0  0
##   0.7  5  4 30 83 72 26  9  1  0
##   0.8 22 28 54 52 37 26  8  2  1
##   0.9 82 39 40 42 25 25 19 16  6
biasx(jj, x_true = 0.8)
## [1] -0.0522
biasr(jj, r_true = 4)
## [1] 0.063
msex(jj, x_true = 0.8)
## [1] 0.0202
mser(jj, r_true = 4)
## [1] 2.977

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  7
##   0.3  0  0  0  0  1  0  0
##   0.4  0  1  0  1  4  0  0
##   0.5  0  0  3  7  5  1  0
##   0.6  0  0  3 13  4  5  0
##   0.7  0  1  9 10  6  0  0
##   0.8  3  2  2  4  4  0  0
##   0.9  6  1  1  0  2  0  1
jj2
##      r
## x      1  2  3  4  5  6
##   0.3  0  0  0  0  1  0
##   0.4  0  0  0  3  2  1
##   0.5  0  2  1  8  6  0
##   0.6  0  0  1 13  5  3
##   0.7  0  4  4 10  1  1
##   0.8  4  3  9  8  0  0
##   0.9  3  2  1  1  3  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.14  1.764  -0.042 0.0218
##  [2,]   -0.31  1.526  -0.034 0.0224
##  [3,]   -0.11  1.386  -0.040 0.0204
##  [4,]   -0.15  1.770  -0.046 0.0182
##  [5,]   -0.26  1.676  -0.034 0.0206
##  [6,]   -0.14  1.384  -0.035 0.0197
##  [7,]   -0.28  1.868  -0.022 0.0170
##  [8,]   -0.21  1.786  -0.014 0.0222
##  [9,]   -0.39  1.994  -0.024 0.0216
mean(fish[,1])
## [1] -0.2211111
mean(fish[,2])
## [1] 1.683778
sd(fish[,1])
## [1] 0.09492687
mean(fish[,2])/sqrt(100)
## [1] 0.1683778
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   8
##   0.3   0   0   1   4   3   2   0   0
##   0.4   0   0   1  10  19   9   0   0
##   0.5   0   0  14  60  55  15   0   0
##   0.6   0   4  40 112  88  20   0   0
##   0.7   7  15  60 113  63  12   1   0
##   0.8  30  40  49  28  17   5   2   0
##   0.9  40  13  13  12  14   2   6   1
jj2
##      r
## x       1   2   3   4   5   6   7   8
##   0.3   0   0   1   3   4   0   0   0
##   0.4   0   0   2  10  23   5   0   0
##   0.5   0   3  20  53  42  14   0   0
##   0.6   0   4  35  96  84  14   0   0
##   0.7   6  19  64 102  72   9   2   0
##   0.8  35  35  43  47  21   5   2   1
##   0.9  35  32  19  21   7   5   3   2
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.133 1.7618 -0.0336 0.01926
##  [2,]  -0.208 1.7688 -0.0212 0.01934
##  [3,]  -0.197 1.7622 -0.0272 0.01952
##  [4,]  -0.202 1.7172 -0.0260 0.01996
##  [5,]  -0.179 1.7554 -0.0226 0.01860
##  [6,]  -0.158 1.6428 -0.0277 0.01909
##  [7,]  -0.148 1.6408 -0.0362 0.02004
##  [8,]  -0.158 1.7108 -0.0281 0.01827
##  [9,]  -0.216 1.7516 -0.0255 0.02035
mean(fish[,1])
## [1] -0.1776667
mean(fish[,2])
## [1] 1.723489
sd(fish[,1])
## [1] 0.02959307
mean(fish[,2])/sqrt(1000)
## [1] 0.0545015

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.00050  0.023700 0.00050 0.0026250
## 0.2 -0.01450  0.001225 0.02200 0.0066775
## 0.3 -0.04750 -0.006050 0.10400 0.0097100
## 0.4 -0.04725 -0.014250 0.25425 0.0120500
## 0.5 -0.11575 -0.014725 0.51075 0.0149475
## 0.6 -0.16400 -0.020250 1.03400 0.0173700
## 0.7 -0.20275 -0.023625 1.75075 0.0191475
## 0.8 -0.00500 -0.054775 2.79500 0.0208225
## 0.9  0.61025 -0.110075 4.61475 0.0271875
matplot(jj,type="b")

cbind(`bias(r)`=jj[,1], `SE(bias)` = jj[,3]/sqrt(4000))
##      bias(r)     SE(bias)
## 0.1 -0.00050 7.905694e-06
## 0.2 -0.01450 3.478505e-04
## 0.3 -0.04750 1.644384e-03
## 0.4 -0.04725 4.020045e-03
## 0.5 -0.11575 8.075667e-03
## 0.6 -0.16400 1.634898e-02
## 0.7 -0.20275 2.768179e-02
## 0.8 -0.00500 4.419283e-02
## 0.9  0.61025 7.296560e-02
cbind(`bias(x)`=jj[,2], `SE(bias)` = jj[,4]/sqrt(4000))
##       bias(x)     SE(bias)
## 0.1  0.023700 4.150489e-05
## 0.2  0.001225 1.055805e-04
## 0.3 -0.006050 1.535286e-04
## 0.4 -0.014250 1.905272e-04
## 0.5 -0.014725 2.363407e-04
## 0.6 -0.020250 2.746438e-04
## 0.7 -0.023625 3.027486e-04
## 0.8 -0.054775 3.292326e-04
## 0.9 -0.110075 4.298721e-04

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.609 -0.0988  0.887 0.02626
## 0.2   0.192 -0.0746  0.980 0.02250
## 0.3  -0.070 -0.0451  1.428 0.02033
## 0.4  -0.255 -0.0291  1.763 0.01899
## 0.5  -0.186 -0.0209  1.926 0.02079
## 0.6   0.057 -0.0224  1.689 0.01870
## 0.7   0.092 -0.0246  1.432 0.02014
## 0.8  -0.062 -0.0433  0.874 0.02045
## 0.9  -0.361 -0.0740  0.511 0.02448
matplot(jj,type="b")

cbind(`bias(r)`=jj[,1], `SE(bias)` = jj[,3]/sqrt(1000))
##     bias(r)   SE(bias)
## 0.1   0.609 0.02804940
## 0.2   0.192 0.03099032
## 0.3  -0.070 0.04515732
## 0.4  -0.255 0.05575096
## 0.5  -0.186 0.06090547
## 0.6   0.057 0.05341087
## 0.7   0.092 0.04528382
## 0.8  -0.062 0.02763831
## 0.9  -0.361 0.01615924
cbind(`bias(x)`=jj[,2], `SE(bias)` = jj[,4]/sqrt(1000))
##     bias(x)     SE(bias)
## 0.1 -0.0988 0.0008304141
## 0.2 -0.0746 0.0007115125
## 0.3 -0.0451 0.0006428910
## 0.4 -0.0291 0.0006005165
## 0.5 -0.0209 0.0006574375
## 0.6 -0.0224 0.0005913459
## 0.7 -0.0246 0.0006368827
## 0.8 -0.0433 0.0006466858
## 0.9 -0.0740 0.0007741256

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.165 -0.0331  1.647 0.02003
## 20   -0.249 -0.0058  1.127 0.01136
## 50   -0.192  0.0090  0.652 0.00620
## 100  -0.123  0.0100  0.301 0.00342
## 200  -0.074  0.0085  0.100 0.00157
matplot(jj,type="b")

cbind(`bias(r)`=jj[,1], `SE(bias)` = jj[,3]/sqrt(1000))
##     bias(r)    SE(bias)
## 10   -0.165 0.052082713
## 20   -0.249 0.035638869
## 50   -0.192 0.020618050
## 100  -0.123 0.009518456
## 200  -0.074 0.003162278
cbind(`bias(x)`=jj[,2], `SE(bias)` = jj[,4]/sqrt(1000))
##     bias(x)     SE(bias)
## 10  -0.0331 6.334042e-04
## 20  -0.0058 3.592347e-04
## 50   0.0090 1.960612e-04
## 100  0.0100 1.081499e-04
## 200  0.0085 4.964776e-05