Manifold horse race model

horse_model.py

import numpy as np
import pymc as pm

# Race times for 6 racers doing 3 races
race_times = np.array([
    [30,20,22],  # Lane 1
    [27,27,28],  # Lane 2
    [25,29,28],  # Lane 3
    [25,23,23],  # Lane 4
    [24,24,25],  # Lane 5
    [20,18,28]   # Lane 6
])

with pm.Model() as model:
    # Priors for racer performance
    mu = pm.Normal('mu', mu=24, sigma=3, shape=6)
    sigma = pm.HalfNormal('sigma', sigma=3, shape=6)
    
    # Reshape mu and sigma to match race_times shape
    mu_reshaped = mu[:, np.newaxis]
    sigma_reshaped = sigma[:, np.newaxis]
    
    # Likelihood of observed race times
    obs = pm.Normal('obs', mu=mu_reshaped, sigma=sigma_reshaped, observed=race_times)
    
    # Posterior distribution of racer performance
    trace = pm.sample(2000, tune=1000)

# Samples from the posterior distribution
samples = trace.posterior['mu'].values

print("Shape of samples:", samples.shape)

# Probability of each racer winning the 4th race
win_probs = np.zeros(6)
for chain in samples:
    for sample in chain:
        win_probs[np.argmin(sample)] += 1

win_probs /= samples.size // 6

print("Likelihood of each lane winning the 4th race:")
for i, prob in enumerate(win_probs):
    print(f"Lane {i+1}: {prob:.4%}")

# Likelihood of each lane winning the 4th race:
# Lane 1: 22.8000%
# Lane 2: 0.2750%
# Lane 3: 1.1375%
# Lane 4: 18.8000%
# Lane 5: 5.3750%
# Lane 6: 51.6125%