Glastocast: predicting the Glastonbury 2027 lineup¶
Live site: glastocast.pages.dev
The question. Which acts get a big slot at Glastonbury 2027 (23 to 27 June)? "Big slot" means anywhere on the Pyramid Stage or the top five on the Other Stage each day, roughly 40 acts a year.
Talking points if I only get five minutes
- The signal is the European summer festival circuit. Acts route one tour through Werchter, Roskilde, Primavera and Glasto, so a booking at one is evidence for the others. Raw effects are confounded by fame, so everything is fame-adjusted.
- The famous fan heuristic, "look for a gap in the tour over Glasto weekend", doesn't work (likelihood ratio about 0.9). A gig in another continent that weekend is the only real red flag.
- My first model double-counted: a prior fitted without festival information, multiplied by odds ratios measured against acts with zero festivals. A forward calibration plot exposed it.
- I once "found" that leave-one-year-out validation flattered my model (13 vs 10 hits a year). Rechecking one change at a time, the scheme barely matters here (11.0 vs 11.3). The gap came from comparing two different models. I'd changed two things at once and blamed the wrong one.
- A hierarchical Bayesian logistic model (partially pooled festival effects, artist effects, random-walk year intercept), fitted with NUTS, beats the best GLM by 2.0 hits a year forward (paired SE 0.45). XGBoost buys nothing: about 550 positives and an additive, monotone signal.
import json, math, re, unicodedata, pickle, warnings, inspect
from collections import defaultdict, Counter
import numpy as np, pandas as pd, matplotlib.pyplot as plt, statsmodels.api as sm
from sklearn.metrics import roc_auc_score
warnings.filterwarnings('ignore')
pd.set_option('display.width', 160); pd.set_option('display.max_columns', 20)
BLUE, ORANGE, AQUA, YELLOW, GREY, RED = '#2a78d6', '#eb6834', '#1baf7a', '#eda100', '#8a877f', '#e34948'
plt.rcParams.update({'figure.dpi': 110, 'axes.spines.top': False, 'axes.spines.right': False, 'axes.grid': True,
'grid.alpha': .25, 'axes.titleweight': 'bold', 'axes.titlesize': 11, 'font.size': 9.5})
def norm(s):
s = unicodedata.normalize('NFKD', s).encode('ascii', 'ignore').decode().lower()
s = re.sub(r'\((band|singer|musician|rapper|group|duo)[^)]*\)', '', s).replace('&', 'and').replace('+', 'and')
s = re.sub(r'^the ', '', s.strip()); return re.sub(r'[^a-z0-9]', '', s)
sig = lambda z: 1 / (1 + np.exp(-z))
RECOMPUTE_SLOW = False # True re-runs the ~12 minute forward MCMC test instead of loading its saved results
1. Data¶
Unit of analysis: one artist in one year, $i=(a,y)$. Label $Y_i = 1$ if they got a big Glasto slot that year.
The candidate pool for year $y$ is everyone who played at least one of ~45 tracked festivals that year. That's the realistic comparison set: acts active on the summer circuit.
Sources
- Glastonbury lineups 2000 to 2025, stage and running order, from Wikipedia.
- Other festivals' lineups from Wikipedia (2000 to 2025, patchy) and setlist.fm (2016 to 2025, much fuller).
- Gig diaries from the setlist.fm API for the weeks around each Glasto weekend.
- Bookmaker odds and 2027 announcements, gathered by hand.
The main trap is lookahead. Final lineups and gig calendars are known after the fact. Anything I use as a feature has to be something you could have known before the Glasto lineup was announced, typically in March. Other festivals announce from November to February, which is why they're usable.
G = json.load(open('data/glasto_top.json'))
heads = pd.DataFrame([r for r in G if r['stage'] == 'Pyramid' and r['rank'] == 0])
heads['artist'] = heads.artist.str.replace(r' \((band|singer|musician)[^)]*\)', '', regex=True)
print('Pyramid headliners, last six editions:')
print(heads.groupby('year').artist.apply(lambda s: ' / '.join(s)).tail(6).to_string())
pool = pd.read_pickle('data/pool3.pkl'); pool = pool[pool.in_pool == 1].reset_index(drop=True)
cfg = json.load(open('data/pool3_cfg.json')); FESTS = cfg['fests']; REGION = cfg['region']
print(f"\npool: {len(pool):,} artist-years, {pool.glasto.sum()} big Glasto slots, {pool.key.nunique():,} artists, {len(FESTS)} festivals")
Pyramid headliners, last six editions: year 2017 Radiohead / Foo Fighters / Ed Sheeran 2019 Stormzy / The Killers / The Cure 2022 Billie Eilish / Paul McCartney / Kendrick Lamar 2023 Arctic Monkeys / Guns N' Roses / Elton John 2024 Dua Lipa / Coldplay / SZA 2025 The 1975 / Neil Young / Olivia Rodrigo pool: 23,731 artist-years, 574 big Glasto slots, 11,649 artists, 46 festivals
# how much setlist.fm adds over Wikipedia, for the festivals that matter most
canon = json.load(open('data/canon.json')); F = json.load(open('data/fest_links_raw.json'))
wiki, sl = defaultdict(set), defaultdict(set)
for f, y, a, _ in F:
m = canon.get(a)
if m and m['mus']: wiki[(f, y)].add(norm(m['c']))
for r in map(json.loads, open('data/fest_setlist.jsonl')):
for nm, mb in r['acts']: sl[(r['fest'], r['year'])].add(norm(nm))
rows = [dict(festival=f, wikipedia=len(wiki.get((f, 2025), ())), setlist_fm=len(sl.get((f, 2025), ())),
both=len(wiki.get((f, 2025), set()) & sl.get((f, 2025), set())))
for f in ['Werchter', 'Roskilde', 'Opener', 'Hurricane', 'MadCool', 'NOSAlive', 'Primavera', 'Pinkpop']]
pd.DataFrame(rows).set_index('festival').T
| festival | Werchter | Roskilde | Opener | Hurricane | MadCool | NOSAlive | Primavera | Pinkpop |
|---|---|---|---|---|---|---|---|---|
| wikipedia | 9 | 5 | 6 | 0 | 0 | 20 | 130 | 25 |
| setlist_fm | 149 | 146 | 70 | 88 | 67 | 52 | 143 | 61 |
| both | 8 | 5 | 6 | 0 | 0 | 19 | 81 | 24 |
2. The festival signal, and why you have to adjust for fame¶
Naively, compare $P(Y=1 \mid \text{played } F)$ with the base rate. But big acts play everything, so every festival looks predictive. The fix is to estimate each festival's effect holding fame and Glasto history fixed:
$$\operatorname{logit} P(Y_i=1) = \alpha_{y} + \beta_f\, f_i + \beta_L L_i + \beta_O O_i + \gamma_F\, b_{iF}$$
- $f_i = \log\big(1 + \sum_{(F',y') : y' < y} 2^{-(y-y')/5}\big)$ is recency-weighted fame: big-festival sets before year $y$, with a 5-year half-life.
- $L_i$ means they played the previous Glasto; $O_i$ means they played an older one.
- $\alpha_y$ are year fixed effects, needed because the pool roughly doubles over the period.
$e^{\gamma_F}$ is the fame-adjusted odds ratio for playing festival $F$ that summer.
BASE = ['fame', 'played_last', 'prev_older']
YD = pd.get_dummies(pool.year.astype(str), prefix='y', drop_first=True, dtype=float)
base_rate = pool.glasto.mean(); res = []
for f in FESTS:
n, hits = int(pool[f].sum()), int(pool.loc[pool[f] == 1, 'glasto'].sum())
if n < 100 or hits < 5: continue
X = sm.add_constant(pd.concat([pool[BASE], YD, pool[[f]]], axis=1))
m = sm.Logit(pool.glasto, X).fit(disp=0)
res.append(dict(festival=f, region=REGION[f], n=n, hits=hits, raw_lift=pool.loc[pool[f] == 1, 'glasto'].mean() / base_rate,
adj_or=np.exp(m.params[f]), lo=np.exp(m.conf_int().loc[f, 0]), hi=np.exp(m.conf_int().loc[f, 1])))
FE = pd.DataFrame(res).sort_values('adj_or').reset_index(drop=True)
fig, ax = plt.subplots(figsize=(7.2, 0.24 * len(FE) + 1))
yy = np.arange(len(FE))
ax.hlines(yy, FE.lo, FE.hi, color=BLUE, lw=1.5, alpha=.5)
ax.scatter(FE.raw_lift, yy, facecolor='none', edgecolor=GREY, s=26, label='raw lift (fame not controlled)', zorder=3)
ax.scatter(FE.adj_or, yy, color=BLUE, s=26, label='fame-adjusted odds ratio, 95% CI', zorder=4)
ax.axvline(1, color='k', lw=.8, ls='--'); ax.set_xscale('log'); ax.set_yticks(yy, FE.festival)
ax.set_xticks([.25, .5, 1, 2, 4, 8], ['0.25', '0.5', '1', '2', '4', '8'])
ax.set_xlabel('effect on chance of a big Glasto slot (log scale)'); ax.legend(loc='lower right', frameon=False)
ax.xaxis.set_minor_formatter(plt.NullFormatter()); ax.set_title('Confounding cuts both ways'); plt.tight_layout(); plt.show()
Fame confounding cuts both ways:
- Festivals that mostly book famous acts look too good raw. Fuji Rock, V Festival and Isle of Wight (grey circle to the right of the blue dot).
- Festivals with deep lineups look too weak raw. Primavera, Roskilde, Flow and Øya book hundreds of small acts, which drags their raw Glasto rate down. Once you compare like with like, Primavera goes from 1.3× to about 2.2×.
What survives adjustment is geography and timing. UK and European June/July festivals sit at 2 to 8×, US festivals at about 2×, and Download (a metal festival a fortnight before Glasto) is the one clear negative.
Bookings across festivals are correlated, since one tour means several festivals. The simplest summary is the number of "circuit" festivals played: UK or European June/July festivals with an adjusted odds ratio of at least 1.5.
SUMMER = {'Primavera','RockAmRing','Pinkpop','Hurricane','Southside','Werchter','Roskilde','Opener','Benicassim','MadCool','NOSAlive','BBK','SuperBock',
'IsleOfWight','Download','Parklife','BST','HydeParkCalling','Wireless','TinthePark','TRNSMT','Latitude','Oxegen','BigWeekend'} # late May to July
CIRC = [r.festival for r in FE.itertuples() if r.adj_or >= 1.5 and r.festival in SUMMER]
pool['n_circ'] = pool[CIRC].sum(axis=1)
g = pool.groupby(pool.n_circ.clip(upper=4)).glasto.agg(['mean', 'size', 'sum'])
fig, ax = plt.subplots(figsize=(5.2, 3))
bars = ax.bar([str(i) if i < 4 else '4+' for i in g.index], g['mean'] * 100, color=BLUE, width=.6)
for b, (p, n) in zip(bars, zip(g['mean'], g['size'])): ax.text(b.get_x() + b.get_width() / 2, p * 100 + 1, f'{p:.0%}\nn={n:,}', ha='center', fontsize=8.5)
ax.set_ylim(0, g['mean'].max() * 125); ax.set_xlabel('circuit festivals played that summer'); ax.set_ylabel('% who got a big Glasto slot')
ax.set_title('Count the festivals'); ax.grid(axis='x', visible=False); plt.tight_layout(); plt.show()
print('circuit festivals:', ', '.join(CIRC))
circuit festivals: SuperBock, HydeParkCalling, Southside, Latitude, Wireless, IsleOfWight, Primavera, MadCool, Benicassim, Hurricane, Parklife, TRNSMT, Pinkpop, NOSAlive, Roskilde, Werchter, Opener, BST, BigWeekend, TinthePark, Oxegen
3. Testing the fan theory: "look for a gap in the tour"¶
Every year, fans look for a hole over Glasto weekend in a big act's tour dates. I pulled setlist.fm gig diaries for the six weeks around Glasto weekend for every big Glasto act from 2010 to 2025. I compared them with the 45 most similar acts each year who didn't play (matched on the festival model's score). That's 980 artist-years.
The lookahead problem. Final calendars include the Glasto set itself, so every Glasto act would trivially show a "gap". I remove the Glasto gig and classify only the other gigs on Friday to Sunday. Then
$$\mathrm{LR}_E = \frac{P(E \mid Y=1)}{P(E \mid Y=0)}$$
for each weekend pattern $E$, with add-one smoothing.
gp = pd.read_pickle('data/gaps.pkl'); core = gp.core_uk + gp.core_eu + gp.core_far
def lr(mask):
a = gp[mask]; pg = (a.glasto.sum() + 1) / (gp.glasto.sum() + 2); png = ((a.glasto == 0).sum() + 1) / ((gp.glasto == 0).sum() + 2)
return dict(n=len(a), played_glasto=round(a.glasto.mean(), 2), LR=round(pg / png, 2))
pd.DataFrame({
'clear gap over the weekend': lr(gp.touring_eu & (core == 0)),
'European gig same weekend': lr((gp.core_eu >= 1) & (gp.core_far == 0)),
'one UK/Irish gig same weekend': lr((gp.core_uk == 1) & (gp.core_far == 0)),
'two+ UK/Irish gigs': lr(gp.core_uk >= 2),
'gig outside Europe that weekend': lr(gp.core_far >= 1),
'no gigs at all nearby': lr((gp.eu_win == 0) & (gp.far_win == 0)),
}).T
| n | played_glasto | LR | |
|---|---|---|---|
| clear gap over the weekend | 416.0 | 0.43 | 0.91 |
| European gig same weekend | 263.0 | 0.42 | 0.88 |
| one UK/Irish gig same weekend | 116.0 | 0.59 | 1.79 |
| two+ UK/Irish gigs | 22.0 | 0.41 | 0.88 |
| gig outside Europe that weekend | 50.0 | 0.14 | 0.22 |
| no gigs at all nearby | 66.0 | 0.55 | 1.46 |
EU = set('GB IE FR DE BE NL LU ES PT IT CH AT DK SE NO FI PL CZ SK HU RO HR SI RS BG GR IS EE LV LT UA TR MT CY BA ME MK AL'.split())
import datetime as dt
cnt = {0: np.zeros(45), 1: np.zeros(45)}; nn = {0: 0, 1: 0}
for r in map(json.loads, open('data/gigs.jsonl')):
fri = dt.date.fromisoformat(r['fri']); nn[r['glasto']] += 1; days = set()
for g_ in r['gigs']:
if g_['city'] in ('Pilton', 'Glastonbury') or 'Worthy' in g_['venue'] or 'Pyramid' in g_['venue'] or g_['cc'] not in EU: continue
o = (dt.datetime.strptime(g_['d'], '%d-%m-%Y').date() - fri).days
if -21 <= o <= 23: days.add(o)
for o in days: cnt[r['glasto']][o + 21] += 1
x = np.arange(-21, 24)
fig, ax = plt.subplots(figsize=(8, 3))
ax.axvspan(-2.5, 2.5, color='#1f6b4a', alpha=.12, lw=0)
ax.plot(x, cnt[1] / nn[1] * 100, color=BLUE, lw=2, label=f'Glasto acts (n={nn[1]})')
ax.plot(x, cnt[0] / nn[0] * 100, color=ORANGE, lw=2, label=f'lookalikes who didn\'t play (n={nn[0]})')
ax.set_xticks(range(-21, 24, 7), ['-21d', '-14d', '-7d', 'Glasto Fri', '+7d', '+14d', '+21d'])
ax.set_ylabel('% with a European gig that day'); ax.legend(frameon=False, loc='upper left')
ax.text(0, ax.get_ylim()[1] * 0.97, 'Glasto Wed to Sun', ha='center', va='top', fontsize=8)
ax.set_title('Everyone gigs at weekends. Glasto acts warm up midweek, then go quiet'); plt.tight_layout(); plt.show()
Verdict. A clear gap is worth about 0.9×, which is no signal at all. Glasto acts are often busy that weekend anyway, doing Dublin or Werchter one night and Worthy Farm the next. A single UK or Irish gig that weekend is actually a good sign (about 1.8×), because acts route a UK run around Glasto. The only real red flag is a gig in another continent (about 0.2×).
4. Model v1, and the bug a calibration plot found¶
v1 was a Bayesian-flavoured update. Start from a prior fitted without festival information, then multiply the odds by a festival-count likelihood ratio:
$$\text{odds}(Y \mid x, n) \;=\; \underbrace{\text{odds}(Y \mid x)}_{\text{prior: bookings unknown}} \times \underbrace{e^{\delta_n}}_{\text{OR of } n \text{ circuit festivals vs } n=0}$$
The bug: $e^{\delta_n}$ came from a joint model where it's the odds ratio relative to acts who played zero circuit festivals. But the prior already averages over everyone, including the many acts with $n \ge 1$. So the prior odds are higher than the $n=0$ odds, and multiplying by $e^{\delta_n}$ counts the festival information twice.
The fix is to fit the conditional directly. For $k = 0, 1, 2, 3$, fit $\operatorname{logit} P(Y \mid x, n \ge k)$ on the rows with $n \ge k$. Conditioning on $n \ge k$ rather than $n = k$ has a second benefit: in October most lineups aren't out yet, and seeing $k$ bookings only tells you the final count is at least $k$.
Both versions get the same per-year level calibration, an intercept shift so predicted totals match recent years, to make the comparison fair.
YEARS = sorted(pool.year.unique()); TEST = [y for y in YEARS if y >= 2007]
pos = pool.groupby('year').glasto.sum()
def shift_to(z, target):
lo, hi = -10, 10
for _ in range(60):
mid = (lo + hi) / 2
lo, hi = (lo, mid) if sig(z + mid).sum() > target else (mid, hi)
return mid
def glm(d, w, cols=BASE):
return sm.GLM(d.glasto, sm.add_constant(d[cols]), family=sm.families.Binomial(), freq_weights=w).fit().params
lin = lambda p, d, cols=BASE: p['const'] + d[cols].values @ p[cols].values
out = []
for y in TEST:
tr, te = pool[pool.year < y], pool[pool.year == y].copy()
w = 0.5 ** ((y - tr.year.values) / 8); target = pos[pos.index < y].tail(5).mean()
nc_tr, nc_te = tr.n_circ.clip(upper=3), te.n_circ.clip(upper=3).values
prior = glm(tr, w); z0 = lin(prior, te); s = shift_to(z0, target)
# v1: prior + log OR(n vs n=0) from a joint model
J = sm.add_constant(pd.concat([tr[BASE], pd.get_dummies(nc_tr, prefix='c', dtype=float).drop(columns='c_0')], axis=1))
jm = sm.GLM(tr.glasto, J, family=sm.families.Binomial(), freq_weights=w).fit().params
z_v1 = z0 + s + np.array([jm.get(f'c_{n}', 0.0) if n else 0.0 for n in nc_te])
# v2: P(Y | x, n >= k) fitted directly
z_v2 = z0.copy()
for k in (1, 2, 3):
m = (nc_tr >= k).values; zk = lin(glm(tr[m], w[m]), te); z_v2 = np.where(nc_te == k, zk, z_v2)
z_v2 = z_v2 + s
out.append(pd.DataFrame(dict(year=y, glasto=te.glasto.values, v1=sig(z_v1), v2=sig(z_v2))))
FW = pd.concat(out)
bins = [0, .01, .02, .05, .1, .2, .35, .5, 1]
fig, ax = plt.subplots(figsize=(5, 4.2))
for col, c, lab in [('v1', RED, 'v1 (double-counted)'), ('v2', BLUE, 'v2 (fixed)')]:
b = FW.groupby(pd.cut(FW[col], bins)).agg(p=(col, 'mean'), o=('glasto', 'mean'), n=('glasto', 'size')).query('n >= 15')
ax.plot(np.sqrt(b.p), np.sqrt(b.o), 'o-', color=c, lw=2, ms=5, label=lab)
t = np.array([0, .01, .05, .1, .25, .5, .75]); ax.plot([0, .9], [0, .9], ls='--', color='k', lw=.8)
ax.set_xticks(np.sqrt(t), [f'{v:.0%}' for v in t]); ax.set_yticks(np.sqrt(t), [f'{v:.0%}' for v in t])
ax.set_xlabel('predicted (square-root scale)'); ax.set_ylabel('actually got a big slot'); ax.legend(frameon=False)
ax.set_title('Forward calibration, 2007 to 2025'); plt.tight_layout(); plt.show()
for col in ('v1', 'v2'):
hits = FW.groupby('year').apply(lambda d: d.nlargest(40, col).glasto.sum()).mean()
print(f'{col}: top-40 hits/yr {hits:.1f} log-loss {-(FW.glasto*np.log(FW[col]) + (1-FW.glasto)*np.log(1-FW[col])).mean():.4f}')
v1: top-40 hits/yr 11.3 log-loss 0.1024 v2: top-40 hits/yr 11.3 log-loss 0.0798
v1 sits well below the diagonal: when it said 40%, those acts played about 10% of the time. But look at the hit rates: identical. More festival bookings still means more likely, so the ordering is right; v1 just overstates by how much. If I'd only tracked top-40 hits, I'd never have known. Calibration plots and log-loss catch bugs that ranking metrics hide.
5. Validation: change one thing at a time¶
My first backtest was leave-one-year-out (fit on all other years, predict the held-out one) and said 13 of 40 right a year. A later forward-chained test (fit only on years before $t$) said about 10. I concluded leave-one-year-out leaked the future. There were two suspects: the scheme itself, and an early "fame" feature, $\log(1 + \#\text{festival sets in all other years})$, which counts later years too.
Same model, changing one thing at a time:
def v2_scores(tr, te, y):
w = 0.5 ** ((y - tr.year.values) / 8); nc_tr, nc_te = tr.n_circ.clip(upper=3), te.n_circ.clip(upper=3).values
z = lin(glm(tr, w), te)
for k in (1, 2, 3):
m = (nc_tr >= k).values; z = np.where(nc_te == k, lin(glm(tr[m], w[m]), te), z)
return z
apps_all = defaultdict(set)
for (f, y), s_ in list(wiki.items()) + list(sl.items()):
for k in s_: apps_all[k].add((f, y))
leaky = pool.copy(); leaky['fame'] = [math.log1p(sum(1 for (_, yy) in apps_all[k] if yy != y)) for k, y in zip(pool.key, pool.year)]
res = {}
for name, d, scheme in [('leaky fame, leave-one-year-out', leaky, 'loyo'), ('clean fame, leave-one-year-out', pool, 'loyo'), ('clean fame, forward-chained', pool, 'fwd')]:
h = []
for y in TEST:
te = d[d.year == y]; tr = d[d.year != y] if scheme == 'loyo' else d[d.year < y]
h.append(te.glasto.values[np.argsort(-v2_scores(tr, te, y))[:40]].sum())
res[name] = round(np.mean(h), 1)
pd.Series(res, name='top-40 hits per year')
leaky fame, leave-one-year-out 11.1 clean fame, leave-one-year-out 11.0 clean fame, forward-chained 11.3 Name: top-40 hits per year, dtype: float64
Neither leaks meaningfully here. All three land within half a hit. Fame changes slowly, and with about 20 years each held-out year is a small slice of the data. The 13 vs 10 gap came from comparing two different models: the early backtest used a flat GLM with all 40 festival flags, the later one the circuit-count model. Section 6 shows the flat GLM really is about 1 hit better when full lineups are known. I'd changed the model and the validation scheme at the same time and blamed the wrong one.
I still use forward-chaining throughout, because it's the only scheme that matches how the model is used: predicting next year from the past.
6. Architecture bake-off: GLM vs XGBoost vs hierarchical¶
Same features, forward-chained from 2007, two scenarios:
- Full lineups known. What you'd have in June.
- Only 30% of bookings announced. Roughly where we are in October. Flat models are trained on randomly masked lineups so they learn that "not announced" isn't the same as "not playing".
(This runs on the Wikipedia-only pool from compare.py, about 15 seconds.)
import subprocess, sys
out = subprocess.run([sys.executable, 'compare.py'], capture_output=True, text=True).stdout # separate process so its globals don't clobber ours
print('\n'.join(out.strip().split('\n')[-9:]))
scenario model hits40 hits_recent auc logloss full lineups A structured GLM (current) 11.27 9.0 0.902 0.0730 full lineups B flat GLM, all 40 flags 12.20 9.6 0.901 0.0716 full lineups C XGBoost, all 40 flags 11.13 9.0 0.901 0.0710 full lineups D hierarchical Bayes (MAP) 11.93 8.6 0.896 0.0722 30% announced A structured GLM (current) 9.00 7.2 0.858 0.0793 30% announced B flat GLM, all 40 flags 8.00 6.4 0.864 0.0801 30% announced C XGBoost, all 40 flags 7.47 7.0 0.843 0.0824 30% announced D hierarchical Bayes (MAP) 8.87 6.4 0.854 0.0808
Everything lands within about 1 hit of everything else, and the year-to-year standard error is about ±0.7. XGBoost buys nothing. There are about 550 positives, and the signal is additive and monotone (more fame, more festivals, more Glasto history), so there are no interactions for trees to find. With partial lineups it's the worst, because trees handle "unknown coded as zero" badly. The structured model's "at least $k$" conditioning handles partial lineups best.
So architecture isn't where the gain is. Pooling information properly is.
7. The hierarchical Bayesian model¶
$$ \begin{aligned} Y_i &\sim \text{Bernoulli}\big(\sigma(\eta_i)\big) \\ \eta_i &= \alpha_{y(i)} + \boldsymbol\beta^\top \mathbf{x}_i + \sum_F \gamma_F\, b_{iF} + u_{a(i)} \\[4pt] \gamma_F &= \mu_{r(F)} + \tau z_F, \quad z_F \sim \mathcal N(0,1) && \text{festival effects, partially pooled by region} \\ \mu_r &\sim \mathcal N(0,1), \quad \tau \sim \text{HalfNormal}(0.7) \\ \alpha_{t} &= \alpha_0 + \sigma_y \textstyle\sum_{s \le t} \epsilon_s, \quad \epsilon_s \sim \mathcal N(0,1) && \text{random-walk year intercept} \\ u_a &= \sigma_a w_a, \quad w_a \sim \mathcal N(0,1) && \text{artist effect} \\ \boldsymbol\beta &\sim \mathcal N(0, 2^2), \quad \alpha_0 \sim \mathcal N(-4, 2^2), \quad \sigma_y \sim \text{HalfNormal}(0.5), \quad \sigma_a \sim \text{HalfNormal}(1) \end{aligned} $$
Why each piece is there:
- Partial pooling. Thin festivals (Open'er, Mad Cool, BST) borrow strength from their region instead of me hand-assigning their effect. Well-measured festivals keep their own estimate.
- Artist effect. Some acts just are Glasto acts, beyond what fame explains. It also soaks up the "metal acts never play Glasto" pattern that a Download flag was standing in for.
- Random-walk year intercept. The pool grows from about 300 to 2,000 acts a year, so the base rate drifts. A random walk means a new year starts from the last one. My first version used independent year effects and predicted 80 Glasto acts a year instead of 28.
- Non-centred parameterisation throughout ($\gamma = \mu + \tau z$ and so on) to avoid funnel geometry when $\tau$ or $\sigma_a$ are small.
Fitted with NUTS in NumPyro (4 chains × 500 draws).
import bayes, numpyro, jax
from numpyro.diagnostics import summary
post, years, arts, mc = bayes.fit(bayes.df, method='nuts', samples=500, warmup=500, chains=4, seed=0)
grouped = mc.get_samples(group_by_chain=True)
s = summary({k: np.asarray(grouped[k]) for k in ('mu', 'tau', 'beta', 'sigma_y', 'sigma_a', 'a0', 'gamma')}, group_by_chain=True)
print('divergences:', int(np.asarray(mc.get_extra_fields()['diverging']).sum()))
pd.DataFrame({k: dict(max_rhat=float(np.max(v['r_hat'])), min_ess=int(np.min(v['n_eff']))) for k, v in s.items()}).T
divergences: 0
| max_rhat | min_ess | |
|---|---|---|
| mu | 1.000825 | 958.0 |
| tau | 1.000746 | 814.0 |
| beta | 1.003861 | 668.0 |
| sigma_y | 1.000841 | 1164.0 |
| sigma_a | 1.006574 | 290.0 |
| a0 | 0.998817 | 1705.0 |
| gamma | 1.001420 | 1963.0 |
g = np.asarray(post['gamma']); order = np.argsort(g.mean(0))
keep = [i for i in order if pool[FESTS[i]].sum() >= 60]
rc = {'EU': BLUE, 'UK': ORANGE, 'US': AQUA, 'AS': YELLOW}
fig, ax = plt.subplots(figsize=(7, 0.22 * len(keep) + 1))
for j, i in enumerate(keep):
lo, hi = np.exp(np.percentile(g[:, i], [5, 95])); c = rc[REGION[FESTS[i]]]
ax.hlines(j, lo, hi, color=c, lw=2, alpha=.55); ax.scatter(np.exp(g[:, i].mean()), j, color=c, s=24, zorder=3)
ax.axvline(1, color='k', lw=.8, ls='--'); ax.set_xscale('log'); ax.set_yticks(range(len(keep)), [FESTS[i] for i in keep])
ax.set_xticks([.5, 1, 2, 4], ['0.5', '1', '2', '4']); ax.set_xlabel('posterior odds multiplier (mean, 90% interval)')
for r, c in rc.items(): ax.scatter([], [], color=c, label={'EU': 'Europe', 'UK': 'UK & Ireland', 'US': 'US', 'AS': 'Asia/Aus'}[r])
ax.xaxis.set_minor_formatter(plt.NullFormatter()); ax.legend(frameon=False, loc='lower right'); ax.set_title('Festival effects, partially pooled by region'); plt.tight_layout(); plt.show()
b = np.asarray(post['beta'])
print(pd.DataFrame({'posterior mean': b.mean(0), 'sd': b.std(0)}, index=['fame (per log-unit)', 'played last Glasto', 'played an older Glasto']).round(2))
print(f"artist effect sd sigma_a = {np.asarray(post['sigma_a']).mean():.2f}, year random-walk step sigma_y = {np.asarray(post['sigma_y']).mean():.2f}, pooling sd tau = {np.asarray(post['tau']).mean():.2f}")
posterior mean sd fame (per log-unit) 1.12 0.09 played last Glasto -0.91 0.30 played an older Glasto 0.48 0.17 artist effect sd sigma_a = 0.66, year random-walk step sigma_y = 0.39, pooling sd tau = 0.40
Reading the coefficients:
- Fame is the biggest single driver.
- Played the last Glasto is negative once the artist effect is in. Given an act's own propensity, they tend to take a year off after playing. In the raw data this looked like no effect, because "Glasto regulars" and "just played" cancel out.
- Played an older Glasto is positive: the comeback effect.
Forward test: Bayesian vs the best GLM¶
Each year from 2007 is refitted from scratch with NUTS on earlier years only, then predicted with a fresh random-walk step for the new year. That's about 12 minutes, so by default this loads the saved run from forward_bayes.py.
if RECOMPUTE_SLOW:
subprocess.run([sys.executable, 'forward_bayes.py'], check=True)
R = pd.read_csv('data/forward_bayes.csv')
d = R.bayes_hits - R.struct_hits
fig, ax = plt.subplots(figsize=(7.5, 3))
ax.plot(R.year.astype(str), R.bayes_hits, 'o-', color=BLUE, lw=2, label=f'hierarchical Bayes (mean {R.bayes_hits.mean():.1f})')
ax.plot(R.year.astype(str), R.struct_hits, 'o-', color=ORANGE, lw=2, label=f'structured GLM (mean {R.struct_hits.mean():.1f})')
ax.set_ylabel('real big-slot acts in top 40'); ax.legend(frameon=False); ax.set_ylim(0, 22)
ax.set_title('Forward test: each year predicted from earlier years only'); plt.tight_layout(); plt.show()
print(f'paired difference: {d.mean():+.2f} hits/yr, sd {d.std():.2f}, se {d.std()/np.sqrt(len(d)):.2f} -> t = {d.mean()/(d.std()/np.sqrt(len(d))):.1f}')
print(f'AUC: Bayes {R.bayes_auc.mean():.3f} vs GLM {R.struct_auc.mean():.3f}')
print(f'predicted Glasto acts per year {R.bayes_sum.mean():.1f} vs actual {R.actual.mean():.1f}, with no hand calibration')
paired difference: +2.00 hits/yr, sd 1.73, se 0.45 -> t = 4.5 AUC: Bayes 0.915 vs GLM 0.909 predicted Glasto acts per year 33.9 vs actual 28.4, with no hand calibration
The Bayesian model is ahead in most years. The average gain is 2.0 hits a year with a paired standard error of about 0.45, and its probability level comes out about right without the hand-tuned calibration the GLM needed.
8. The 2027 board¶
The live site currently runs the structured GLM:
- prior times the at-least-$k$ festival update;
- bookie headline prices as a floor, with the margin removed (minimum 1.6×, because the press only reports the top of each book);
- tour-weekend likelihood ratios;
- bootstrap ranges (refitting on resampled years, alternative bookies, and the slot-count assumption);
- a multinomial stage model.
Swapping in the Bayesian model is next: posterior predictive with unannounced bookings marginalised out.
html = open('site/index.html').read()
D = json.loads(re.search(r'const D=(\{.*?\});\n', html, re.S).group(1))
B = pd.DataFrame([dict(act=a['name'], p=a['p'], lo=a['lo'], hi=a['hi'],
bookies=a['bookie']['odds'] if a['bookie'] else '',
headline_if_booked=a['stage']['head'], other_stage_if_booked=a['stage']['oth'],
evidence='; '.join(e.get('short') or e['label'] for e in a['evidence'])) for a in D['acts'][:20]])
B.index = range(1, 21)
B.style.format({'p': '{:.0%}', 'lo': '{:.0%}', 'hi': '{:.0%}', 'headline_if_booked': '{:.0%}', 'other_stage_if_booked': '{:.0%}'})
| act | p | lo | hi | bookies | headline_if_booked | other_stage_if_booked | evidence | |
|---|---|---|---|---|---|---|---|---|
| 1 | Tame Impala | 57% | 46% | 59% | 9% | 52% | tour gap; Roskilde, Open'er, Rock Werchter, NOS Alive | |
| 2 | Sam Fender | 45% | 38% | 46% | 5/4 | 62% | 25% | |
| 3 | Harry Styles | 36% | 19% | 37% | 4/6 | 90% | 2% | tour gap |
| 4 | Robbie Williams | 30% | 28% | 32% | EVS | 90% | 3% | tour gap |
| 5 | Fontaines D.C. | 29% | 25% | 32% | 4/1 | 44% | 43% | 2 UK/IE gigs that weekend |
| 6 | Florence and the Machine | 28% | 25% | 31% | 4/1 | 44% | 38% | |
| 7 | Muse | 28% | 23% | 31% | 8/1 | 69% | 17% | tour gap; Hurricane / Southside |
| 8 | Arctic Monkeys | 27% | 24% | 30% | 6/1 | 58% | 36% | |
| 9 | Foals | 25% | 21% | 29% | 12% | 71% | ||
| 10 | Little Simz | 24% | 21% | 27% | 6/1 | 37% | 27% | tour gap |
| 11 | Fred Again | 22% | 20% | 26% | 6/1 | 40% | 37% | |
| 12 | Becky Hill | 21% | 17% | 24% | 6% | 63% | ||
| 13 | Foo Fighters | 20% | 16% | 23% | 55% | 23% | ||
| 14 | Two Door Cinema Club | 20% | 16% | 23% | 9% | 52% | ||
| 15 | The Killers | 20% | 16% | 23% | 31% | 42% | ||
| 16 | Caroline Polachek | 20% | 16% | 24% | 13% | 46% | Primavera Sound, Roskilde | |
| 17 | Raye | 19% | 15% | 24% | 5/1 | 54% | 22% | EU gig that weekend; Open'er |
| 18 | Royal Blood | 19% | 15% | 22% | 12% | 37% | ||
| 19 | Ed Sheeran | 19% | 17% | 24% | 5/1 | 56% | 17% | |
| 20 | Charli XCX | 18% | 16% | 25% | 6/1 | 49% | 43% | no Euro dates |
9. Limitations and what I'd do next¶
- Coverage. About 30% of Glasto's big acts play none of the tracked festivals that year, so this kind of model can never see them coming. Reunions, legends and surprise sets mostly sit there.
- Bookies and evidence overlap. If the bookies have already priced in a festival booking, multiplying by its likelihood ratio counts it twice. The Bayesian version should treat the bookie price as an observation with its own noise, not as a floor.
- Partial lineups in October. The proper treatment is a booking model, $P(b_{iF} = 1 \mid x_i, \text{announced})$, and marginalising over unannounced festivals. The at-least-$k$ trick is a decent approximation.
- Tour evidence was estimated on a matched set (Glasto acts plus lookalikes), not the full pool, and assumed independent of festival bookings given $Y$.
- Genre from Wikidata is noisy (Robbie Williams is tagged "electronic music"), which limits the stage model (63% right vs a 53% baseline).