Modelling time series where model output feeds into data

Questions
1

I’m trying to model a rating system for a competitive game that comprises of 2 parts: a match prediction network and a rating update network, each having several parameters.

Each team starts off with the same rating.
Iterating through the matches in ascending chronological order, a logistic regression predicts the match outcome based on a single predictor: the difference in the 2 team’s ratings, so that there is an intercept (b0) and a coefficient (b1) to be identified.
The predicted outcome from the logistic regression and the actual outcome are then used with a gain parameter (alpha) to detail how much rating should be assigned to the winner (and consequently removed from the loser).
Each team has their rating accordingly updated.
After several rounds of competition the ratings have steadied out.

i.e. for game i in 1…N:
team1 = games[i][‘team1’]
team2 = games[i][‘team2’]
logit(Pr(team1win)) = b0 + b1 * (rating[team1] - rating[team2])
rating_update = alpha * (outcome[i] - Pr(team1win))
ratings[team1] += rating_update
ratings[team2] -= rating_update

I’ve framed it as a general optimisation problem by using the total accuracy over the dataset as the fitness value and fitting the 3 parameters using a genetic algorithm works well, but I’d be interested in treating the parameters as random variables using PyMC.
However, I can’t work out how to parameterise it since the data inputs (ratings and Pr(team1win)) aren’t fixed in advance and are dependent upon the output of the 2 models at each timestep.

2

Instead of looking at individual rounds I’d just write this down as one logistic regression.

with pm.Model():
    raw = pm.Normal('rating_raw', sd=1, shape=n_teams)
    rating = pm.Deterministic('rating', rating_raw * rating_sd)
    pm.Deterministic('elo', 1400 + 400 * rating / np.log(10))
    pm.Bernoulli('winner', logit_p=rating[team1] - rating[team2], observed=team1_won)

where team1 and team2 are arrays with n_games entries, that each contain the index of the corresponding team in game i. team1_won is 1 if team1[i] won game i and 0 otherwise.

2 Likes
3

Thanks for the reply!
I’ve implemented this model and adapted it slightly, as in the third line both ‘rating_raw’ and ‘rating_sd’ were undefined. What is the purpose of this line actually, I don’t understand the difference in the standard deviation of rating_raw and multiplying rating_raw by a (different) sd?

with pm.Model() as mod_logistic:
    rating_raw = pm.Normal('rating_raw', sd=1, shape=n_teams)
    rating = pm.Deterministic('rating', rating_raw * 2)
    pm.Deterministic('elo', 1400 + 400 * rating / np.log(10))
    pm.Bernoulli('winner', logit_p=rating[team1] - rating[team2], observed=team1_won)

On my dataset (I’ll try and upload a simulated test one soon) it has worked and produced sensible ratings. However, these ratings are averages over the timeframe covered by the training data. The reason I implemented this as a time series in my Genetic Algorithm method was so that I can track each how each team’s rating changes over time.

image
image.jpg888×438 217 KB

4

I forgot to define those vars…

with pm.Model() as mod_logistic:
    rating_sd = pm.HalfNormal('rating_sd', sd=2)
    rating_raw = pm.Normal('rating_raw', sd=1, shape=n_teams)
    rating = pm.Deterministic('rating', rating_sd * rating_raw)
    pm.Deterministic('elo', 1400 + 400 * rating / np.log(10))
    pm.Bernoulli('winner', logit_p=rating[team1] - rating[team2], observed=team1_won)

You could also write that as

with pm.Model() as mod_logistic:
    rating_sd = pm.HalfNormal('rating_sd', sd=2)
    rating = pm.Normal('rating', sd=rating_sd, shape=n_teams)
    pm.Deterministic('elo', 1400 + 400 * rating / np.log(10))
    pm.Bernoulli('winner', logit_p=rating[team1] - rating[team2], observed=team1_won)

The first one is a centered, the second a non-centered parametrization. They give the same results, but the sampler usually likes one of them much more than the other.
A nice explanation is here.

About the time dependence. That is I think a bit more tricky, but with a bit of experimentation it should be doable.

First, I guess you have to pick a definition for “time”. Is it the number of games played? Actual time? The log of the game number (people usually improve faster at the beginning)? I’d start with a simple linear change of rating over time (whichever). While not particularly realistic, it should be pretty simple to implement and gives you something to compare other models to. Something like that:

with pm.Model() as mod_logistic:
    rating_mu_sd = pm.HalfNormal('rating_sd', sd=2)
    rating_mu = pm.Normal('rating', sd=rating_mu_sd, shape=n_teams)

    sd = pm.HalfNormal('rating_time_sd', sd=2)
    raw = pm.Normal('rating_time', sd=1, shape=n_teams)
    rating_time = pm.Deterministic('rating_time', sd * raw)

    rating_time1 = rating_time[team1] * team1_log_gametime
    rating_time2 = rating_time[team2] * team2_log_gametime
    theta = (rating[team1] + rating_time1) - (rating[team2] + rating_time2)

    max_rating = rating + max_log_gametime
    pm.Deterministic('max_elo', 1400 + 400 * max_rating / np.log(10))
    pm.Bernoulli('winner', logit_p=theta, observed=team1_won)

Where team1_log_gametime is an array/dataframe of shape (n_games) that contains in at index i the log of how many games team1 played up to that particular game. max_log_gametime contains the log of the number of games each team has played up to now. You might want to standardize all log_gametime values so that mean(max_log_gametime) == 0 and std(max_log_gametime) == 1.

It’s probably also a good idea to think a bit about the mean of rating_time. Does it even make sense to have negative values here? Somehow we are measuring the performance of teams not absolute but in relation to other teams, so maybe?

If you are lucky the linear log_gametime works mostly. If you want to look at non-linear changes over time you could also try GPs, but with more than a few teams that probably takes some work to get right.

1 Like
5

That’s really useful, I hadn’t thought about including time in that manner.
The only issue now is that I don’t follow what max_log_gametime is and am getting an input mismatch where it should be added to rating_mu, where the latter has nteams on dimension 0 and max_log_gametime has ngames dimension 0. This is because I thought it should be an ngames vector with each entry i being max(team1_log_gametime[i], team2_log_gametime[i]), but I can’t work out from your description or from playing around what it should be.

with pm.Model() as mod_time:
    rating_mu_sd = pm.HalfNormal('rating_mu_sd', sd=2)
    rating_mu = pm.Normal('rating_mu', sd=rating_mu_sd, shape=n_teams)

    time_sd = pm.HalfNormal('time_sd', sd=2)
    time_raw = pm.Normal('time_raw', sd=1, shape=n_teams)
    rating_time = pm.Deterministic('rating_time', time_sd * time_raw)

    rating_time1 = rating_time[team1] * team1_log_gametime
    rating_time2 = rating_time[team2] * team2_log_gametime
    theta = (rating_mu[team1] + rating_time1) - (rating_mu[team2] + rating_time2)

    max_rating = rating_mu + max_log_gametime
    pm.Deterministic('max_elo', 1400 + 400 * max_rating / np.log(10))
    pm.Bernoulli('winner', logit_p=theta, observed=team1_won)
6

Does that help?

df = pd.DataFrame({
    'team1': ['blueberry', 'bilberry', 'strawberry', 'blueberry'],
    'team2': ['bilberry', 'strawberry', 'blueberry', 'strawberry'],
    'team1_won': [False, True, False, True],
    'time': ['2019-07-01', '2019-07-02', '2019-07-03', '2019-07-04'],
    'team1_gamenum': [1, 2, 2, 3], # How many games did team1 already play?
    'team2_gamenum': [1, 1, 2, 3],
}).rename_axis(index='game')
df['time'] = pd.to_datetime(df['time'])
# How many games were played by 
df['team1_log_gamenum'] = np.log(df['team1_gamenum'])
df['team2_log_gamenum'] = np.log(df['team2_gamenum'])

max_log_gametime = np.log(df.team1.value_counts() + df.team2.value_counts())

max_log_gametime.index = max_log_gametime.index.astype('category')
df = df.astype({
    'team1': max_log_gametime.index.dtype,
    'team2': max_log_gametime.index.dtype})
n_teams = len(max_log_gametime)

with pm.Model() as model:
    rating_mu_sd = pm.HalfNormal('rating_mu_sd', sd=2)
    rating_mu = pm.Normal('rating_mu', sd=rating_mu_sd, shape=n_teams)

    sd = pm.HalfNormal('rating_time_sd', sd=2)
    raw = pm.Normal('rating_time_raw', sd=1, shape=n_teams)
    rating_time = pm.Deterministic('rating_time', sd * raw)

    rating_time1 = rating_time[df.team1.cat.codes] * df['team1_log_gamenum']
    rating_time2 = rating_time[df.team2.cat.codes] * df['team2_log_gamenum']
    theta = (
        (rating_mu[df.team1.cat.codes] + rating_time1)
        - (rating_mu[df.team2.cat.codes] + rating_time2))

    max_rating = rating_mu + max_log_gametime.values
    pm.Deterministic('max_elo', 1400 + 400 * max_rating / np.log(10))
    pm.Bernoulli('winner', logit_p=theta, observed=df.team1_won)
    
    trace = pm.sample()
2 Likes
7

Yes it does indeed, thank you so much for your help!

1 Like