Calculating loglikelihood for permutations / race outcomes

Modeling
1

Hi all,

I am reading a paper about a hierarchical model that uses race data.

For each race, j, they observe the finishing position of driver i. From these finishing positions, they estimate the racing ability of every driver. They start with a model where the racing ability of driver i is constant across each race:

image

and

image

In order to create a finishing order of the drivers from the racing abilities, they begin with the last driver. They choose a driver to finish in last place with probability proportional to

image

Next, they choose the second-to-last driver from the remaining drivers again with probability proportional to lambda_i,j This process is repeated until drivers i,1 and i,2 remain. The probability that driver i,1 finished second is:

image

In the step before, there were three drivers left in the equation. In this case, the third place driver i,1 is chosen with probability:

image

The probability that driver i,2 becomes third is:

image

Basically, in the model drivers drop out one by one until the last two drivers are left. Intuitively, each driver remains in the race for an exponential length of time with rate (inverse mean) parameter, lambda_i,j. The longer a driver remains in the race, the better the racing ability of this driver.

First question: how exactly is the the last-place driver selected? We know that the probability that a variable is the minimum of a number of exponential random variables is:

image

In the paper the authors say the following:

“We choose a driver to finish in last place with probability proportional to image”

Let’s say there are 20 drivers. Why not compare the probabilities of each driver getting last as follows:

image

and see which driver has the highest likelihood to finish last. Its not clear to me how they classify one the 20 drivers as the last place driver

Next, let image be defined so that image means that driver i finishes in k-th place in race j. If we have J races and I drivers then we have the following likelihood function:

image

I am trying to reproduce this in Stan-code:

data {
int<lower=1> n_drivers;
int<lower=1> n_races;
int<lower=1, upper=n_drivers> driver_id[n_races];
matrix[n_drivers, n_races] race_results // matrix with the finishing positions of each driver for all the races
}

parameters {
real<lower=0> tau_race_skills;
vector[n_drivers - 1] raw_race_skills;
}

transformed parameters {
vector[n_drivers] race_skills;
vector[n_drivers] lambda_race_skills

for (driver in 1:(n_drivers - 1)) {
race_skills[driver] = raw_race_skills[driver];
}

race_skills[n_drivers] = -sum(raw_race_skills);
lambda_race_skills = exp(-race_skills);
}

model {
// priors
raw_race_skills ~ normal(0, tau_race_skills);

// likelihood

}

How do I write the loglikelihood function in my stan code? In addition, how do I make sure that it follows the method where each driver drops out one-by-one?

Help and/or tips are much appreciated!!

3 
data { 
int<lower=1> n_drivers; 
int<lower=1> n_races;
int<lower=1, upper=n_drivers> driver_id[n_races];
matrix[n_races, n_drivers] race_results
}

parameters {
real<lower=0> tau_race_skills;
vector[n_drivers - 1] raw_race_skills;
}

transformed parameters {
vector[n_drivers] race_skills;
vector[n_drivers] lambda_race_skills

for (driver in 1:(n_drivers - 1)) {
    race_skills[driver] = raw_race_skills[driver];
}

race_skills[n_drivers] = -sum(raw_race_skills);    
lambda_race_skills = exp(-race_skills);
}

model {  
// priors 
raw_race_skills ~ normal(0, tau_race_skills);

// likelihood
for (race in 1:n_races) {
    for (driver in 1:n_drivers) {
        // determine kth place
        finish_position = n_drivers - driver + 1
        // find driver_id of driver i that finished kth place in race j 
        driver = race_results[race, finish_position]
        // find driver_id of drivers that didn't drop out yet in race j
        other_drivers = race_results[race, 1:finish_position]
        lambda_race_skills_other_drivers = 0
        for (other_driver in 1:other_drivers) {
            lambda_race_skills_other_drivers += lambda_race_skills[other_driver]
        }
        // calculate prob that driver i finished kth place in race j 
        target += lambda_race_skills[driver] / lambda_race_skills_other_drivers
    }  
}

Is the above Stan-interpretation of the paper methodology correct? Please advice! Help is much appreciated!

4

Any time you have a custom probability you want to use, you can use the target += syntax to add the log probability directly to the log density that Stan samples.

So this shows how you can do that for a normal distribution: 7.3 Increment log density | Stan Reference Manual

For a bernoulli with probability p, you could replace:

y ~ bernoulli(p); // assume y is integer 0 or 1

with:

if(y == 0) {
  target += log(1 - p);
} else {
  target += log(p);
}

(though the second will be slower and isn’t vectorized)

Here is a paper that might be related to the problem you are looking at: http://www.glicko.net/research/multicompetitor.pdf .

I think this video might be related to the problem you’re working with as well: Statistical Rethinking Winter 2019 Lecture 14 - YouTube