Recommendation, Ranking, and Personalization Questions

The idea

Both methods start from the same table of who interacted with what, and both score an unseen item by borrowing from neighbours. User-based collaborative filtering finds users whose history resembles yours and recommends what they liked. Item-based collaborative filtering finds items resembling the ones you already interacted with, where "resembling" means the same people interacted with both. The difference is which axis the similarity is computed on, and that decides cost, how dense the data behind each similarity is, and how well the method scales.

Mechanics and cost

For a table with $U$ users and $I$ items, the similarity table to build has one entry per pair:

• User-based: about $U(U-1)/2$ user pairs.
• Item-based: about $I(I-1)/2$ item pairs.

Sparse implementations do not compare all pairs. They count co-occurrences: each user with $m_u$ interactions contributes about $m_u^2$ item pairs to the item-item table, and each item with $n_i$ interactions contributes about $n_i^2$ user pairs to the user-user table. The user-user cost is therefore dominated by popular items (an item seen by a million users contributes about $10^{12}$ user pairs), while the item-item cost is dominated by heavy users, who can be capped or sampled.

A tiny trace. Three users and four items (1 means the user interacted with the item). Cosine similarity treats each row (or column) as a list of numbers and measures the angle between two of them: the dot product divided by the product of the lengths, which for 0/1 data is 0 when nothing is shared and 1 when two histories are identical. The program scores user 0's unseen items 2 and 3 both ways, using all neighbours since the table is tiny (the real systems below keep only the k most similar, "neighbours"):

import numpy as np

R = np.array([[1, 1, 0, 0],      # user 0 has seen items 0 and 1
              [1, 1, 1, 0],      # user 1
              [0, 1, 1, 1]], float)

def cosine(M):
    n = np.linalg.norm(M, axis=1)
    return (M @ M.T) / np.outer(n, n)

su = cosine(R)[0]                # similarity of every user to user 0
user_scores = {j: round(float(sum(su[u] * R[u, j] for u in (1, 2))), 3) for j in (2, 3)}
si = cosine(R.T)                 # item-item similarity
item_scores = {j: round(float(sum(si[i, j] for i in (0, 1))), 3) for j in (2, 3)}
print("user-user similarity to user 0 (users 1, 2):", np.round(su[1:], 3))
print("user-based scores:", user_scores)
print("item-item similarity of items 0 and 1 to item 2:", np.round([si[0, 2], si[1, 2]], 3))
print("item-based scores:", item_scores)
m = R.sum(1)
print("item pairs each user contributes (m*(m-1)/2):", (m * (m - 1) / 2).astype(int))

Output:

user-user similarity to user 0 (users 1, 2): [0.816 0.408]
user-based scores: {2: 1.225, 3: 0.408}
item-item similarity of items 0 and 1 to item 2: [0.5   0.816]
item-based scores: {2: 1.316, 3: 0.577}
item pairs each user contributes (m*(m-1)/2): [1 3 3]

Reading it: user-based asks which other users resemble user 0 (user 1 at 0.816 and user 2 at 0.408) and adds up what they interacted with, so item 2 scores $0.816 + 0.408 = 1.225$ and item 3 only $0.408$. Item-based asks which items resemble the ones user 0 already has, so item 2 scores $0.5 + 0.816 = 1.316$ and item 3 scores $0 + 0.577 = 0.577$. Both put item 2 first here; they differ in which table they needed. The last line shows the co-occurrence count: a user with $m$ interactions adds $m(m-1)/2$ item pairs, which grows with $m^2$, so one very heavy user can dominate the item-item cost, as described above.

At 10M users and 1M items, the last program below prints the pair counts: 5.000e+13 user pairs against 5.000e+11 item pairs, a ratio of 100 (users outnumber items 10 to 1, and the pair count grows with the square, so $10^2 = 100$), which is the first reason item-based is the usual production choice when users outnumber items.

Sensitivity to sparsity

The quality of a similarity depends on how many interactions overlap. A cosine similarity (defined above) computed from a zero overlap is just zero, and one computed from a single shared interaction is noisy. Which side is denser depends on the data shape, not on a rule: with many more users than items, each item vector (a column of the table) usually holds many interactions while each user vector holds few, which favours item-based. A long-tail catalog reverses that for the tail: a rarely seen item has a very short vector and its similarities are weak.

The program makes this concrete with 4000 users, 300 items and 12 interactions per user, at two levels of popularity skew (each item gets a popularity weight $e^{\sigma z}$ with $z$ a standard normal draw, a lognormal weight: with $\sigma = 0.5$ the weights are similar, with $\sigma = 2.0$ a few items get huge weights and most get tiny ones). It counts pairs that share nothing and runs a leave-one-out test (hide one interaction per tested user and see whether it is recovered): score unseen items using each user's or item's 50 nearest neighbours (the 50 most similar users or items; all other similarities are set to zero), and check whether the hidden item lands in the top 10 (hit@10). A popularity baseline is included, since a neighbourhood method must beat it to be worth running.

The code builds the 0/1 table X (users in six taste clusters, each pulled toward its own items and weighted by popularity), hides one interaction per tested user, turns similarities into scores with topk_rows (keep only each row's 50 largest similarities) and one matrix product, then ranks the unseen items. A last block repeats the comparison on seeds 1 to 20 without the per-run detail and prints the mean hit@10 of each method, the mean difference between item-based and user-based, and how many seeds each one wins.

import numpy as np

def cos(M):
    n = np.linalg.norm(M, axis=1) + 1e-9
    S = (M @ M.T) / np.outer(n, n); np.fill_diagonal(S, 0); return S

def topk_rows(S, k):
    idx = np.argpartition(-S, k, axis=1)[:, :k]
    out = np.zeros_like(S)
    np.put_along_axis(out, idx, np.take_along_axis(S, idx, axis=1), axis=1)
    return out

def run(popularity_sigma, seed=0, n_users=4000, n_items=300, per_user=12, k=50, verbose=True):
    rng = np.random.default_rng(seed)
    cluster_of_item = rng.integers(0, 6, n_items)
    pop = np.exp(popularity_sigma * rng.normal(size=n_items))      # lognormal popularity skew
    X = np.zeros((n_users, n_items), dtype=bool)
    for u in range(n_users):
        c = rng.integers(0, 6)
        w = pop * np.where(cluster_of_item == c, 6.0, 1.0)
        X[u, rng.choice(n_items, per_user, replace=False, p=w / w.sum())] = True
    Xf = X.astype(np.float32)
    uo, io = Xf @ Xf.T, Xf.T @ Xf
    iu, ii = np.triu_indices(n_users, 1), np.triu_indices(n_items, 1)
    user_empty, item_empty = np.mean(uo[iu] == 0), np.mean(io[ii] == 0)
    if verbose:
        print(f"popularity sigma {popularity_sigma}: items with under 5 interactions: {np.mean(X.sum(0) < 5):.3f}")
        print(f"  user pairs sharing no item: {user_empty:.3f}   item pairs sharing no user: {item_empty:.3f}")
    test_users = rng.choice(n_users, 600, replace=False)
    hidden = np.array([rng.choice(np.where(X[u])[0]) for u in test_users])
    Xtr = X.copy(); Xtr[test_users, hidden] = False
    Xt = Xtr.astype(np.float32)
    ub = topk_rows(cos(Xt)[test_users], k) @ Xt
    ib = Xt[test_users] @ topk_rows(cos(Xt.T), k)
    popscore = np.tile(Xt.sum(0), (len(test_users), 1))
    hits = {}
    for name, sc in [("popularity", popscore), ("user-based", ub), ("item-based", ib)]:
        sc = np.where(Xtr[test_users], -np.inf, sc)
        top = np.argsort(-sc, axis=1)[:, :10]
        hits[name] = np.array([hidden[j] in top[j] for j in range(len(test_users))])
        if verbose:
            print(f"  {name} hit@10: {hits[name].mean():.3f}")
    if verbose:
        print(f"  users where exactly one of user-based / item-based hits: {int((hits['user-based'] != hits['item-based']).sum())} of 600")
    return user_empty, item_empty, {name: h.mean() for name, h in hits.items()}

print("4000 users x 300 items, 12 interactions per user, 50 neighbours")
run(0.5); run(2.0)

print("\nseeds 1 to 20, mean over seeds")
for sigma in (0.5, 2.0):
    out = [run(sigma, seed=sd, verbose=False) for sd in range(1, 21)]
    diff = np.array([h["item-based"] - h["user-based"] for _, _, h in out])
    ahead, behind = int((diff > 1e-9).sum()), int((diff < -1e-9).sum())
    print(f"sigma {sigma}: user pairs sharing no item {np.mean([o[0] for o in out]):.3f}, "
          f"item pairs sharing no user {np.mean([o[1] for o in out]):.3f}")
    print("  hit@10 " + ", ".join(f"{n} {np.mean([o[2][n] for o in out]):.3f}" for n in ("popularity", "user-based", "item-based")))
    print(f"  item-based minus user-based: {100 * diff.mean():+.1f} points; item-based ahead on {ahead} seeds, "
          f"user-based ahead on {behind}, tied on {20 - ahead - behind}")

U, I = 10_000_000, 1_000_000
print(f"pairs to compare at {U:,} users: {U*(U-1)//2:.3e}; at {I:,} items: {I*(I-1)//2:.3e}")
print(f"ratio of the two pair counts: {(U*(U-1)//2)/(I*(I-1)//2):.0f}")

Output:

4000 users x 300 items, 12 interactions per user, 50 neighbours
popularity sigma 0.5: items with under 5 interactions: 0.000
  user pairs sharing no item: 0.538   item pairs sharing no user: 0.062
  popularity hit@10: 0.117
  user-based hit@10: 0.182
  item-based hit@10: 0.182
  users where exactly one of user-based / item-based hits: 58 of 600
popularity sigma 2.0: items with under 5 interactions: 0.170
  user pairs sharing no item: 0.006   item pairs sharing no user: 0.644
  popularity hit@10: 0.475
  user-based hit@10: 0.483
  item-based hit@10: 0.492
  users where exactly one of user-based / item-based hits: 59 of 600

seeds 1 to 20, mean over seeds
sigma 0.5: user pairs sharing no item 0.548, item pairs sharing no user 0.059
  hit@10 popularity 0.093, user-based 0.161, item-based 0.187
  item-based minus user-based: +2.6 points; item-based ahead on 18 seeds, user-based ahead on 0, tied on 2
sigma 2.0: user pairs sharing no item 0.018, item pairs sharing no user 0.609
  hit@10 popularity 0.459, user-based 0.506, item-based 0.495
  item-based minus user-based: -1.0 points; item-based ahead on 4 seeds, user-based ahead on 15, tied on 1
pairs to compare at 10,000,000 users: 5.000e+13; at 1,000,000 items: 5.000e+11
ratio of the two pair counts: 100

Reading it: with mild skew (sigma 0.5) only 0.062 of item pairs share no user but 0.538 of user pairs share no item, so item vectors are the dense side. With heavy skew (sigma 2.0) it reverses: 0.170 of items have under 5 interactions and 0.644 of item pairs share no user, while only 0.006 of user pairs share no item, because every user touches the popular head. So "item-based is less sensitive to sparsity" is true for a catalog without a heavy tail and false for one with it.

Accuracy follows the dense side, by small margins. On the first run the two tie at sigma 0.5 (0.182 each) and are 0.9 points apart at sigma 2.0 (0.483 against 0.492), but one seed is not the picture. Over seeds 1 to 20 the last block of output shows item-based ahead at sigma 0.5 by 2.6 points on average (ahead on 18 seeds, behind on none), and user-based ahead at sigma 2.0 by 1.0 point (ahead on 15 seeds, against 4 for item-based). The two methods are different rather than interchangeable (on the first run they disagree on 58 and 59 of 600 users), yet close, because with 12 interactions per user and 50 neighbours both recover the six taste clusters easily, so accuracy is capped by the data more than by the axis. At sigma 2.0 the popularity baseline reaches 0.475 on the first run (0.459 on average over seeds 1 to 20), within 5 points of either method, which is a warning that heavy skew lets popularity do most of the work. So the case for item-based rests mainly on cost, staleness, precomputed serving and explanations, plus a small accuracy edge when items are the dense side; the honest answer to "which is more accurate" is to test on your own table. The data is synthetic; the numbers show the mechanism, not a ranking of the methods.

When to prefer each

• Item-based: many more users than items, a fairly stable catalog, a need for fast serving from precomputed neighbour lists, and a wish to say why an item was suggested. It is the default for large consumer catalogs.
• User-based: few users and many items (a niche community or enterprise setting), taste that shifts fast so that "people who just did what you did" matters, or when user-user closeness is itself the product, as in social features.
• Beyond both: matrix factorization or two-tower embeddings (models that learn a short vector for each user and each item, so a pair with no shared history can still be scored) replace neighbourhoods when the table gets large and sparse, because they generalize across items that share no users.

Common mistakes

• Using raw overlap without shrinkage: two items with one shared user and cosine 1.0 look perfect. Shrink similarity toward zero when the overlap is small, for example multiply it by $n/(n+\lambda)$ where $n$ is the number of shared users and $\lambda$ a constant you tune (with $\lambda = 10$, one shared user keeps only 1/11 of the similarity).
• Ignoring popularity: popular items are similar to everything. Normalize by item popularity or use a lift-style measure (how much more often two items occur together than their individual popularities would predict).
• Recomputing user neighbours in batch only, so a user's afternoon activity changes nothing until tomorrow.
• Claiming one method is more accurate without testing it on your own table: the output above shows the answer depends on its shape.

Running the code

Save each program above to a file and run it with python3. Python 3.12 and numpy 2.1.3 were used, and seeds are fixed inside the code.

The idea in plain words

A K-armed Bernoulli bandit is a row of K options (arms), for example K article thumbnails. Pulling arm $i$ pays 1 (a click) with an unknown probability $p_i$ and 0 otherwise. The goal is to collect as many clicks as possible while you are still learning which arm is best. That creates a tension: playing the arm that looks best now (exploitation) versus trying others to learn (exploration).

Thompson sampling resolves it with one idea: keep a belief about each $p_i$, draw one random plausible value from each belief, and play the arm whose draw is largest. Arms you are unsure about produce widely spread draws, so they still win sometimes; arms you are confident are bad almost never win. Exploration shrinks automatically as evidence accumulates. Epsilon-greedy instead explores a fixed fraction $\varepsilon$ of the time (uniformly random arm) and otherwise plays the arm with the best average so far, so it keeps wasting a fixed share of traffic forever.

Mechanism

For Bernoulli rewards the natural belief is a Beta distribution, a distribution over a probability between 0 and 1. Start every arm at Beta(1, 1) (uniform, "no idea"). After arm $i$ has produced $s_i$ successes and $f_i$ failures, its posterior (the belief after seeing data) is exactly

Each round: (1) draw $\theta_i \sim \text{Beta}(1 + s_i, 1 + f_i)$ for every arm, (2) play $\arg\max_i \theta_i$, (3) observe the reward $r \in \{0, 1\}$ and add it to that arm's success or failure count. Everything is O(K) per round (the work per round grows in proportion to the number of arms, since each arm needs one draw) and stores two numbers per arm. The notation $\arg\max_i \theta_i$ means "the index $i$ of the arm whose draw is largest". The update works because the Beta family is conjugate to Bernoulli data: a model whose posterior has the same form as its prior, so each outcome just adds 1 to one of the two numbers.

Regret is the yardstick. Pseudo-regret after $T$ rounds is $\sum_{t=1}^{T} (p^* - p_{a_t})$, where $p^* = \max_i p_i$ and $a_t$ is the arm played in round $t$: the clicks you gave up compared with always playing the best arm. It is computed from the true rates (known only in a simulation), not from the noisy rewards, which removes luck from the comparison.

A by-hand example: an arm that has produced 3 clicks in 10 tries has posterior Beta(4, 8), with mean 4/12 = 0.333 and a wide spread (standard deviation 0.131), so its draws land anywhere from about 0.15 to 0.6. The same observed rate after 300 clicks in 1000 tries gives Beta(301, 701), with mean 0.300 and standard deviation 0.014, so its draws sit within a few hundredths of 0.3. The first program lines print both, with four draws from each. The wide belief is what lets an under-tested arm win a draw now and then.

Implementation and comparison

Setup: five arms with rates 0.05, 0.10, 0.12, 0.15, 0.20 (the gaps are small, as click rates usually are), 5000 rounds, 100 independent seeds per algorithm, epsilon-greedy with $\varepsilon = 0.1$ that pulls each arm once before it starts. A seed fixes the random number generator, so each of the 100 runs is a repeatable, independent replay of the whole 5000-round game, and averaging over 100 runs smooths out luck. The two algorithms use different seed ranges, so their runs are independent.

import numpy as np
from scipy.stats import beta as beta_dist

def run_ts(p, T, rng):
    K = len(p)
    a = np.ones(K); b = np.ones(K)          # Beta(1,1) prior per arm
    best = p.max(); regret = np.zeros(T); r_sum = 0.0
    for t in range(T):
        arm = int(np.argmax(rng.beta(a, b)))  # sample one plausible rate per arm, play the largest
        r = float(rng.random() < p[arm])
        a[arm] += r; b[arm] += 1 - r
        r_sum += best - p[arm]
        regret[t] = r_sum
    return regret

def run_eps(p, T, rng, eps=0.1):
    K = len(p)
    n = np.zeros(K); s = np.zeros(K)
    best = p.max(); regret = np.zeros(T); r_sum = 0.0
    exploit_loss = np.zeros(T)                   # regret paid on steps where the exploit branch ran
    for t in range(T):
        if t < K:
            arm = t                              # pull every arm once first
        elif rng.random() < eps:
            arm = int(rng.integers(K))           # explore: uniform random arm
        else:
            arm = int(np.argmax(s / n))          # exploit: best empirical mean
            exploit_loss[t] = best - p[arm]
        r = float(rng.random() < p[arm])
        n[arm] += 1; s[arm] += r
        r_sum += best - p[arm]
        regret[t] = r_sum
    return regret, exploit_loss

p = np.array([0.05, 0.10, 0.12, 0.15, 0.20])
T, seeds = 5000, 100
ts = np.array([run_ts(p, T, np.random.default_rng(s)) for s in range(seeds)])
runs = [run_eps(p, T, np.random.default_rng(1000 + s)) for s in range(seeds)]
eg = np.array([x[0] for x in runs]); exploit_loss = np.array([x[1] for x in runs])
# What a belief looks like: 3 clicks in 10 tries versus 300 in 1000, both from a Beta(1, 1) start.
for s_, f_ in [(3, 7), (300, 700)]:
    d = beta_dist(1 + s_, 1 + f_)
    print(f"{s_} clicks in {s_ + f_}: Beta({1 + s_}, {1 + f_}) mean={d.mean():.3f} sd={d.std():.3f} "
          f"draws={np.round(d.rvs(4, random_state=0), 3).tolist()}")
print("arms", p.tolist(), "seeds", seeds)
print("horizon  TS_mean  eps0.1_mean  TS_wins_of_100")
for h in (100, 500, 1000, 2500, 5000):
    print(f"{h:7d}  {ts[:, h-1].mean():7.1f}  {eg[:, h-1].mean():11.1f}  {(ts[:, h-1] < eg[:, h-1]).sum():14d}")
# eps-greedy keeps exploring: regret per step late in the run
print("slope per step, steps 1000-2500: TS %.4f  eps %.4f" % (
    (ts[:, 2499] - ts[:, 999]).mean() / 1500, (eg[:, 2499] - eg[:, 999]).mean() / 1500))
print("slope per step, steps 2500-5000: TS %.4f  eps %.4f" % (
    (ts[:, -1] - ts[:, 2499]).mean() / 2500, (eg[:, -1] - eg[:, 2499]).mean() / 2500))
# Where the late eps-greedy slope comes from: random picks versus exploit steps that picked a non-best arm.
gap = p.max() - p
print("late eps-greedy, per step: random picks 0.1 x mean gap = %.4f; exploit steps that missed the best arm = %.4f" % (
    0.1 * gap.mean(), exploit_loss[:, 2500:].mean()))

Output:

3 clicks in 10: Beta(4, 8) mean=0.333 sd=0.131 draws=[0.481, 0.618, 0.282, 0.399]
300 clicks in 1000: Beta(301, 701) mean=0.300 sd=0.014 draws=[0.318, 0.331, 0.296, 0.309]
arms [0.05, 0.1, 0.12, 0.15, 0.2] seeds 100
horizon  TS_mean  eps0.1_mean  TS_wins_of_100
    100      6.2          8.3              65
    500     20.1         23.1              55
   1000     28.6         35.0              56
   2500     40.5         59.4              60
   5000     49.3         91.0              77
slope per step, steps 1000-2500: TS 0.0080  eps 0.0163
slope per step, steps 2500-5000: TS 0.0035  eps 0.0126
late eps-greedy, per step: random picks 0.1 x mean gap = 0.0076; exploit steps that missed the best arm = 0.0049

Reading the result

• Early (100 and 500 rounds) the two are close: regret 6.2 versus 8.3 at 100 rounds and 20.1 versus 23.1 at 500, and Thompson sampling has the lower regret in 65 and 55 of the 100 seeds. With gaps this small, no algorithm can separate the arms from a few hundred pulls, so a single run is nearly a coin flip. This is why a claim that one bandit beats another needs many seeds and several horizons, not one run.
• By 5000 rounds the gap is clear: 49.3 versus 91.0 mean regret, with Thompson sampling ahead in 77 of 100 seeds.
• The slope explains why. Over rounds 2500 to 5000 Thompson sampling gives up about 0.0035 clicks per round and epsilon-greedy about 0.0126. Regret that grows linearly in $T$ means a constant loss per round that never goes away; logarithmic growth means the loss per round keeps shrinking, so each doubling of $T$ adds only a fixed amount. Thompson sampling's slope falls from 0.0080 (rounds 1000 to 2500) to 0.0035, which is the slowing that logarithmic growth predicts (the logarithmic rate itself is a known result for Bernoulli arms; this table shows the slowing, it does not prove the rate). Epsilon-greedy's slope falls only from 0.0163 to 0.0126 and does not go to zero.
• The late epsilon-greedy slope has two parts, both printed. First, a fixed $\varepsilon = 0.1$ sends 10% of traffic to a uniformly random arm forever; a random arm loses on average the mean gap to the best arm, 0.076, so this costs about $0.1 \times 0.076 = 0.0076$ per round. Second, the other 90% of rounds are exploit steps, and with the best arm only 0.05 ahead of the runner-up and 0.08 ahead of the third arm, the arm with the best average is sometimes not the best arm, which costs about 0.0049 per round in rounds 2500 to 5000. The two parts add to the measured 0.0126. A decaying $\varepsilon$ removes the first part and the second shrinks as the averages sharpen; Thompson sampling has no separate knob for either.

Trade-offs and pitfalls

• Epsilon-greedy is easier to explain and has one knob, but the knob is a permanent tax. A decaying $\varepsilon_t$ (for example proportional to $1/t$) fixes the linear growth but needs a schedule that depends on the unknown gaps.
• Thompson sampling needs a posterior. Beta-Bernoulli is closed-form; other reward types need another conjugate model or an approximation. Its randomness also makes single runs noisy, so report averages over seeds.
• Delayed or batched feedback: if clicks arrive hours later, you keep sampling from a stale posterior (a belief that does not yet include recent outcomes, so it is more confident in old information than it should be). Thompson sampling copes better than greedy rules because it keeps randomizing, but update in batches and sample once per request.
• Non-stationary rates (the best arm changes over time): both algorithms need forgetting (discounted counts or a sliding window); otherwise old data dominates.
• A bug to avoid: using the observed reward instead of the true rate when computing regret in a simulation makes the curves noisy and the comparison harder to read.

The idea in plain words

A play, a skip or a thumbs-up is behaviour, not a statement of taste. A play can be started by autoplay; a skip can mean "not now"; a thumbs-up is rare and comes from a self-selected minority (people who chose to give feedback, so they are not a random sample of users). Turning these into labels means deciding, for each event, three things: which direction it points (positive, negative, or unknown), how much to trust it (a weight), and what to leave out (ambiguous events dropped rather than guessed). Raw "played = positive" labels are mostly noise; a layered label with weights is far cleaner.

Label design, signal by signal

Normalize dwell by item length (fraction consumed, not minutes), or long items look more liked than short ones.

Hard negatives, weighting and label bias

• Hard negatives are evidenced rejections: early skips and thumbs-down, where the user saw the item and turned it down. This is the only meaning of the term used below, and the simulation tests this kind. A different source of negatives is model-mined negatives, items the current model scores high that the user ignored; they contain many false negatives (items the user would have liked but did not get to), so if you use them, give them a low weight.
• Weighting. Each example carries a weight $w_i$ in the loss $\sum_i w_i\,\ell(y_i, f(x_i))$. Explicit feedback gets the largest weight, autoplay-driven signals are halved, weak negatives get the smallest.
• Propensity-weighted loss. Signals are observed only for items that were shown, and shown in particular positions. To correct, divide by the estimated chance that the signal could have been observed, $\hat p_i$ (the propensity), clipped so no weight explodes:

where $\ell$ is the per-example loss and $p_{\min}$ is a floor such as 0.1. A propensity is estimated, for example, from a small slice of traffic whose order was randomized: the play rate at each position divided by the play rate at position 1. The by-hand example below works it; the simulation in the next section does not apply this weighting.

• Label bias to watch. Label bias is any systematic reason a label differs from true taste. The common sources are autoplay (the system chose the play), position (top slots get plays regardless of taste), length (long items have lower completion), who bothers to give thumbs (a minority, skewed to extremes), and the old recommender (users only play what it showed).
• Sparse clicks plus a large unlabeled log, with weekly updates. The unlabeled log is mostly "shown, nothing happened". A first version needs two things from this list. First, treat the unlabeled log as weak negatives with a small weight, which is the idea of positive-unlabeled learning: learning when only some positives are labelled and the rest of the data is a mix of positives and negatives, so unlabeled does not mean negative. Second, retrain weekly on a rolling window (only the most recent period of data, so the model follows current taste), warm-started from the previous model (starting from its weights, not from scratch), with a replay of older positives (a sample of older positive examples mixed into each retrain) so rare tastes are not forgotten. An optional add-on when labelled positives are too sparse is to pretrain on the unlabeled log: self-supervised pretraining uses labels that come from the data itself, such as predicting which item a user played next, and sampled-softmax pretraining scores the true item against a random sample of other items instead of the whole catalogue.

Worked example: propensity weights by hand, then how much cleaner the labels are

Propensity weights by hand. From a randomized slice, plays per 1,000 impressions are 300 at position 1, 150 at position 3 and 30 at position 10, so the propensities relative to position 1 are 1.00, 0.50 and 0.10. Each event's base weight from the label table is divided by its propensity (floored at $p_{\min}$). The program prints the result for two floors.

import numpy as np

# Estimate the propensity of a position from a small slice of traffic whose order was randomized:
# plays per 1,000 impressions at each position (illustrative numbers), divided by the rate at position 1.
plays_per_1000 = {1: 300, 3: 150, 10: 30}
p_hat = {pos: n / plays_per_1000[1] for pos, n in plays_per_1000.items()}

# (event, base weight w_i from the label table, position where it was shown)
events = [("user-chosen play", 1.0, 1), ("autoplay completion", 0.5, 3), ("user-chosen play", 1.0, 10), ("thumbs-up", 3.0, 10)]
print("event, position, p_hat, base weight, final weight with p_min 0.1, final weight with p_min 0.25")
for name, w, pos in events:
    p = p_hat[pos]
    print(f"{name}, {pos}, {p:.2f}, {w:.1f}, {w / max(p, 0.1):.1f}, {w / max(p, 0.25):.1f}")

Output:

event, position, p_hat, base weight, final weight with p_min 0.1, final weight with p_min 0.25
user-chosen play, 1, 1.00, 1.0, 1.0, 1.0
autoplay completion, 3, 0.50, 0.5, 1.0, 1.0
user-chosen play, 10, 0.10, 1.0, 10.0, 4.0
thumbs-up, 10, 0.10, 3.0, 30.0, 12.0

A user-chosen play at position 10 counts ten times as much as the same play at position 1 (10.0 against 1.0), because only one tenth as many impressions there lead to a play, and the thumbs-up at position 10 jumps from 3.0 to 30.0. That is a heavy weight for one event, and the floor $p_{\min}$ is the control: raising it from 0.1 to 0.25 cuts those weights to 4.0 and 12.0 at the price of under-correcting the deepest positions. Choose it by checking that the weighted loss helps on a randomized slice rather than only adding variance.

Label construction on a simulation. A simulation of 40,000 impressions. A hidden flag says whether the user would like the item (about 30% would). Autoplay starts the item for 40% of impressions, which makes about three quarters of all plays autoplay (the rest are user-chosen). A user-chosen play is five times as likely for a liked item (0.50 against 0.10); consumed fraction is higher for liked items; thumbs are given on 4% of plays and mostly match taste. The code builds labels as in the table and measures precision (the share of positives that the user really likes) and the share of liked items hiding among the negatives. The hidden flag exists only because this is a simulation; in production it is unobservable.

import numpy as np
import pandas as pd

rng = np.random.default_rng(0)
N = 40000
pref = rng.random(N) < 0.30                                    # latent truth: would the user like this item
autoplay = rng.random(N) < 0.40                                # playback started by the system, not by a user choice
started = autoplay | (rng.random(N) < np.where(pref, 0.50, 0.10))
frac = np.where(pref, rng.beta(5, 2, N), rng.beta(1.5, 4, N))  # fraction of the item consumed if started
frac = np.where(started, frac, 0.0)
thumb = started & (rng.random(N) < 0.04)
thumb_up = thumb & (rng.random(N) < np.where(pref, 0.9, 0.15))
thumb_dn = thumb & ~thumb_up

df = pd.DataFrame(dict(pref=pref, autoplay=autoplay, started=started, frac=frac, up=thumb_up, down=thumb_dn))

def build_labels(d):
    """label 1/0 and a weight; rows that stay NaN are dropped from training instead of guessed."""
    label = pd.Series(np.nan, index=d.index)
    w = pd.Series(0.0, index=d.index)
    done = d.started & (d.frac >= 0.7) & ~d.down
    label[done] = 1; w[done] = np.where(d.autoplay[done], 0.5, 1.0)       # completion counts less when the system chose it
    skip = d.started & (d.frac < 0.1) & ~d.up
    label[skip] = 0; w[skip] = np.where(d.autoplay[skip], 0.5, 1.0)        # early skip: negative, a hard-negative candidate
    label[d.up] = 1; w[d.up] = 3.0                                          # explicit feedback overrides, with the largest weight
    label[d.down] = 0; w[d.down] = 3.0
    weak = ~d.started                                                       # shown, never played: weak negative
    label[weak] = 0; w[weak] = 0.1
    return label, w

label, w = build_labels(df)
naive = df.started
def report(name, pos_mask, neg_mask, weights=None):
    wt = np.ones(len(df)) if weights is None else weights
    tp = (wt * (pos_mask & df.pref)).sum(); pp = (wt * pos_mask).sum()
    fn_in_neg = (wt * (neg_mask & df.pref)).sum() / max((wt * neg_mask).sum(), 1e-9)
    print(f"{name:34s} positives={int(pos_mask.sum()):6d}  precision={tp / pp:.3f}  liked-share-among-negatives={fn_in_neg:.3f}")
print("share of items the user truly likes:", round(df.pref.mean(), 3))
report("naive: played = positive", naive, ~naive)
pos_mask, neg_mask = (label == 1), (label == 0)
report("constructed labels (unweighted)", pos_mask, neg_mask)
report("constructed labels (weighted)", pos_mask, neg_mask, w.values)
hard = (label == 0) & (w >= 0.5)
print("hard negatives (skip or thumbs-down):", int(hard.sum()), " liked share among them:", round(float(df.pref[hard].mean()), 3))
weak_neg = (label == 0) & (w == 0.1)
print("weak negatives (never started):", int(weak_neg.sum()), " liked share among them:", round(float(df.pref[weak_neg].mean()), 3))
print("rows dropped as ambiguous:", int(label.isna().sum()))

Output:

share of items the user truly likes: 0.295
naive: played = positive           positives= 21343  precision=0.387  liked-share-among-negatives=0.190
constructed labels (unweighted)    positives=  5260  precision=0.940  liked-share-among-negatives=0.168
constructed labels (weighted)      positives=  5260  precision=0.925  liked-share-among-negatives=0.103
hard negatives (skip or thumbs-down): 2691  liked share among them: 0.014
weak negatives (never started): 18657  liked share among them: 0.19
rows dropped as ambiguous: 13392

Reading it:

• Treating every play as a positive gives 21,343 positives, of which only 0.387 are items the user likes (autoplay plays the system chose, and plays that were abandoned).
• The constructed labels give 5,260 positives with precision 0.940 (0.925 when weights are applied). That is far fewer positives, but each is far more reliable.
• Hard negatives (skips and thumbs-down) number 2,691 and only 0.014 of them are items the user likes, so they are safe to use at full weight.
• The weak negatives, shown but never started, number 18,657 and 0.190 of them ARE items the user likes. That is why they carry a weight of 0.1 and not 1.0.
• The liked share among all negatives is 0.190 for the naive labels, 0.168 for the constructed labels unweighted and 0.103 once the weights apply. The unweighted gain is small because the 18,657 weak negatives make up most of the negatives; weighting them at 0.1 is what cuts the share of liked items hiding among the negatives from 0.190 to 0.103.
• Skips look safe at full weight here (0.014 liked) partly because of how the simulator is built: a liked item's consumed fraction is drawn so that it is almost never below 10% (about 0.005% of started liked items). In real logs a user may skip an item they like because of mood, context or a duplicate, so the liked share among early skips is not zero, and the cut-off and weight should be checked on held-out data.
• 13,392 ambiguous rows are dropped rather than guessed.

Common mistakes

• Using play or click as a positive with no weighting for autoplay or position.
• Treating "no interaction" as a confirmed negative at full weight.
• Not normalizing dwell time by item length.
• Choosing the completion and skip cut-offs by feel. Tune them against the downstream metric on held-out data, and re-check them when the product changes autoplay behaviour.

Direct answer

Treat each customer arrival as a decision: the context is what you know about the arrival and the stock, the action is a price, and the reward is the margin earned (price $p$ minus unit cost $c$, if the item sells). Inventory turns this into a bandit with a budget: every sale uses up a unit that cannot be sold later. The standard fix is a shadow price $\lambda$ (also called a Lagrange multiplier or the dual variable of the stock constraint: the price, in money, of using up one unit of stock): choose the price that maximises $(p-c-\lambda)\cdot \Pr(\text{sale}\mid x,p)$, and adjust $\lambda$ up when stock is selling faster than the plan and down when slower. Validate offline with logs that contain randomised prices (inverse propensity scoring, doubly robust estimates) plus a demand-model simulator for the inventory path, then roll out behind price bounds, a holdout (a group kept on the current pricing for comparison), and a stock-pace kill switch (an automatic switch back to safe pricing when sales stray from plan).

The design, piece by piece

State and context $x$. Customer segment or intent signals, product features, time left in the selling window, remaining inventory, competitor price if you have it, day of week. Inventory and time-left are part of the state because the same price is right early and wrong late. Keep protected attributes out of the price decision.

Action $a$. A price from a discrete grid (for example steps of a few percent around a reference price). A grid keeps the logs easy to evaluate; see the continuous options below.

Reward. Realised margin $(p-c)\cdot \mathbb{1}[\text{sale}]$ for the arrival, with $c$ the unit cost. Optimise a longer horizon quantity (total margin until the end of the window) through $\lambda$, not by adding a hand-made penalty to each sale.

Constraints.

1. Inventory: total sales over the horizon at most the stock $I$ (the pace target is $I/T$ sales per arrival).
2. Price floor and ceiling, and a maximum price change per step, so prices are not erratic.
3. Margin floor (never below cost plus a minimum).
4. Consistency and fairness rules: no price differences driven by protected attributes, and no more than a set spread between customers at the same moment.
5. Exploration limits: only explore inside a band around the current price.

Learner. A demand model for $q(x,p)=\Pr(\text{sale}\mid x,p)$ with the shape built in: a logistic curve that decreases in price, so evidence at one price informs its neighbours. Thompson sampling draws the demand parameters from their posterior (the distribution of plausible values given the sales so far) and prices as if the draw were true, so uncertain prices get tried; the dual variable is updated by $\lambda \leftarrow \max(0,\ \lambda + \eta(\text{sales}_t - I/T))$. This is the bandit-with-knapsacks idea: a single resource budget priced by a multiplier.

Worked example from a run

The program takes the unit cost as 0, so the reward is the price $p$ when a sale happens (with a cost $c$, replace $p$ by $p - c$ everywhere, so the rule becomes maximise $(p - c - \lambda)\, q$). The demand curve is $q(s,p) = 1 / (1 + e^{-(a_s - 0.2p + 1)})$ for segment $s$, where $a_s$ is 1.0 for segment 0 and 3.0 for segment 1 (the segment's baseline appetite), 0.2 is how fast demand falls with price, and the +1 is a shared offset. This is a known stand-in for a curve you would fit. There are two equal-sized customer segments, 11 prices from 10 to 30, 1000 arrivals and 120 units in stock. The program finds, by bisection, the shadow price at which expected sales cross the pace $I/T = 0.120$ per arrival. Because the price grid is coarse, expected sales jump in steps as $\lambda$ moves, so the crossing is a jump and the program keeps both of its ends: the high end (sales just under the pace) and the low end (sales just over it). It compares always taking the myopic best price ($\lambda=0$) with a policy at each end and with a mixture of the two that hits the pace exactly, and then checks whether logs with uniformly random prices recover the shadow-priced policy's value by inverse propensity scoring (IPS: each logged row is weighted by the target policy's probability of the logged price divided by the logging policy's probability of that price; with continuous prices that denominator is the logging density).

import numpy as np

prices = np.arange(10, 31, 2.0)                     # action set: price grid
seg_a = np.array([1.0, 3.0])                         # two customer segments (the context), 50% each
b = 0.2
T, INV = 1000, 120                                   # periods (one arrival each) and units in stock

def q(s, p):  # purchase probability of segment s at price p
    return 1 / (1 + np.exp(-(seg_a[s] - b * p + 1.0)))
def best_price(s, lam):                              # maximise (p - lam) * q: lam is the price of a unit of stock
    return prices[np.argmax((prices - lam) * q(s, prices))]
def sales_rate(lam): return np.mean([q(s, best_price(s, lam)) for s in (0, 1)])

lo, hi = 0.0, 40.0
for _ in range(40):                                  # find lam so that expected sales per period = INV / T
    mid = (lo + hi) / 2
    lo, hi = (mid, hi) if sales_rate(mid) > INV / T else (lo, mid)
lam, lam_lo = hi, lo                                 # the two ends straddle a jump in expected sales (the price grid is coarse)
w_lo = (INV / T - sales_rate(lam)) / (sales_rate(lam_lo) - sales_rate(lam))   # share of arrivals priced at lam_lo to hit the plan exactly
print(f"target sales rate {INV/T:.3f}; shadow price lam={lam:.2f}; sales rate at lam {sales_rate(lam):.3f}; just below it {sales_rate(lam_lo):.3f}; at lam=0 {sales_rate(0):.3f}")

def simulate(lam_used, seed, lam_alt=None, w_alt=0.0):
    rng = np.random.default_rng(seed); stock = INV; rev = 0.0; sold_out_at = T
    for t in range(T):
        if stock == 0:
            sold_out_at = min(sold_out_at, t); break
        s = rng.integers(2)
        l = lam_alt if (lam_alt is not None and rng.random() < w_alt) else lam_used
        p = best_price(s, l)
        if rng.random() < q(s, p): stock -= 1; rev += p
    return rev, sold_out_at

print("\nrevenue and stock-out period, mean of 200 seeds")
for name, l, alt, wa in (("greedy (lam=0)", 0.0, None, 0.0), ("shadow, high side", lam, None, 0.0),
                         ("shadow, low side", lam_lo, None, 0.0), ("shadow, mixed", lam, lam_lo, w_lo)):
    res = np.array([simulate(l, s, alt, wa) for s in range(200)])
    print(f"{name:18s} revenue {res[:,0].mean():7.1f}   period stock ran out {res[:,1].mean():6.1f} (1000 = never)")

# offline check: log with uniformly random prices, estimate the lagrangian policy's value per period by IPS
def true_value(l): return np.mean([best_price(s, l) * q(s, best_price(s, l)) for s in (0, 1)])
def ips(seed, n=20000):
    rng = np.random.default_rng(seed)
    s = rng.integers(2, size=n); pa = rng.integers(len(prices), size=n); p = prices[pa]
    buyd = rng.random(n) < q(s, p); r = p * buyd
    match = np.array([best_price(si, lam) for si in s]) == p
    w = match * len(prices)                          # target prob 1 on the chosen price, logging prob 1/len(prices)
    return (w * r).mean(), (w * r).sum() / w.sum(), (w * buyd).mean()   # revenue IPS, self-normalised, sales-rate IPS
vals = np.array([ips(sd) for sd in range(20)])
print(f"\nlagrangian policy, 20 logs of 20000 uniformly random prices each")
print(f"revenue per period: true {true_value(lam):.3f}, IPS mean {vals[:,0].mean():.3f} (sd across logs {vals[:,0].std():.3f}), self-normalised {vals[:,1].mean():.3f}")
print(f"sales per period: true {sales_rate(lam):.3f}, IPS mean {vals[:,2].mean():.3f} (sd {vals[:,2].std():.3f})")

Output:

target sales rate 0.120; shadow price lam=20.71; sales rate at lam 0.104; just below it 0.135; at lam=0 0.595

revenue and stock-out period, mean of 200 seeds
greedy (lam=0)     revenue  1615.2   period stock ran out  201.1 (1000 = never)
shadow, high side  revenue  2843.8   period stock ran out  998.9 (1000 = never)
shadow, low side   revenue  3107.4   period stock ran out  887.0 (1000 = never)
shadow, mixed      revenue  3090.4   period stock ran out  965.4 (1000 = never)

lagrangian policy, 20 logs of 20000 uniformly random prices each
revenue per period: true 2.861, IPS mean 2.834 (sd across logs 0.195), self-normalised 2.812
sales per period: true 0.104, IPS mean 0.103 (sd 0.007)

• The myopic policy sells at about 0.595 per arrival and runs out around arrival 201, so roughly 800 later arrivals find an empty shelf; the shadow-priced policies keep sales near the plan and earn 2843.8 to 3107.4 against 1615.2 on average over 200 seeds, depending on which end of the jump described next they use. The gap comes from this demand curve and stock level, and is not a general number.
• The price grid is coarse, so expected sales jump in steps as $\lambda$ moves, and the bisection ends on both sides of one jump: at 20.71 expected sales are 0.104 per arrival, under the plan of 0.120, and a hair lower they are 0.135, over it. The two ends are different policies. The high end leaves about 16 of the 120 units unsold and earns 2843.8; the low end sells nearly all the stock (it runs out at period 887 on average) and earns 3107.4; randomising between them so that expected sales equal 0.120 earns 3090.4 and runs out at period 965.4. So the revenue of "the" shadow-price policy depends on which end of the jump it is built from, and a comparison against it should name the end.
• IPS recovers the high-end policy's value (true 2.861 revenue per period, IPS mean 2.834) and its sales rate (true 0.104, IPS 0.103), but a single log of 20,000 rows has a standard deviation of 0.195 around that value because only one price in eleven matches the policy on each row. The constraint metric (sales pace) must be evaluated offline too, not only revenue.

Seeing the shadow price work, and the online loop

The next program prints, for each segment, the price chosen and the chance of a sale at $\lambda = 0$ and at the shadow price 20.71, then runs the loop from the design with demand unknown: each arrival, Thompson sampling draws a plausible purchase probability for every price from a Beta posterior (a table of 22 segment-price cells stands in for the smooth logistic model, so the program needs no fitting routine), the price maximising $(p - \lambda)\, q$ is shown, the cell is updated with the outcome, and $\lambda \leftarrow \max(0, \lambda + \eta(\text{sale} - I/T))$ with $\eta = 2$.

import numpy as np

prices = np.arange(10, 31, 2.0)
seg_a = np.array([1.0, 3.0])
b = 0.2
T, INV = 1000, 120

def q(s, p):
    return 1 / (1 + np.exp(-(seg_a[s] - b * p + 1.0)))

def best_price(s, lam):
    return prices[np.argmax((prices - lam) * q(s, prices))]

# A by-hand case: what each segment is offered and sells at lam = 0 and at the shadow price.
for lam in (0.0, 20.71):
    print(f"lam = {lam:5.2f}")
    for s in (0, 1):
        p = best_price(s, lam)
        print(f"  segment {s}: price {p:4.1f}  P(sale) {q(s, p):.3f}  (price - lam) * P(sale) = {(p - lam) * q(s, p):.3f}")
    print(f"  expected sales per arrival: {np.mean([q(s, best_price(s, lam)) for s in (0, 1)]):.3f}")

# The online loop from the text, with demand unknown: Thompson sampling on a Beta posterior for every
# (segment, price) cell, and the dual update lam <- max(0, lam + eta * (sale - INV / T)).
def online(seed, eta=2.0):
    rng = np.random.default_rng(seed)
    wins = np.ones((2, len(prices))); losses = np.ones((2, len(prices)))      # Beta(1, 1) priors
    lam, stock, rev, sold_out = 0.0, INV, 0.0, T
    lams = []
    for t in range(T):
        if stock == 0:
            sold_out = t
            break
        s = rng.integers(2)
        draw = rng.beta(wins[s], losses[s])                                  # one plausible demand curve
        k = int(np.argmax((prices - lam) * draw))
        sale = rng.random() < q(s, prices[k])
        wins[s, k] += sale; losses[s, k] += 1 - sale
        if sale:
            stock -= 1; rev += prices[k]
        lam = max(0.0, lam + eta * (float(sale) - INV / T))
        lams.append(lam)
    return rev, sold_out, float(np.mean(lams[len(lams) // 2:]))

res = np.array([online(sd) for sd in range(200)])
print("\nonline learner (eta=2), mean of 200 seeds")
print(f"revenue {res[:, 0].mean():.1f}, period stock ran out {res[:, 1].mean():.1f} (1000 = never), mean lam over the second half of the run {res[:, 2].mean():.2f}")
print(f"seeds where the stock ran out before period 1000: {(res[:, 1] < T).sum()} of 200; earliest {res[:, 1].min():.0f}, latest {res[:, 1].max():.0f}")

Output:

lam =  0.00
  segment 0: price 10.0  P(sale) 0.500  (price - lam) * P(sale) = 5.000
  segment 1: price 16.0  P(sale) 0.690  (price - lam) * P(sale) = 11.040
  expected sales per arrival: 0.595
lam = 20.71
  segment 0: price 26.0  P(sale) 0.039  (price - lam) * P(sale) = 0.207
  segment 1: price 28.0  P(sale) 0.168  (price - lam) * P(sale) = 1.225
  expected sales per arrival: 0.104

online learner (eta=2), mean of 200 seeds
revenue 3013.6, period stock ran out 913.7 (1000 = never), mean lam over the second half of the run 21.18
seeds where the stock ran out before period 1000: 200 of 200; earliest 750, latest 952

Reading it:

• At $\lambda = 0$ segment 1 is offered 16 and buys with probability 0.690, so it earns 11.040 per arrival in the objective; at $\lambda = 20.71$ the same segment is offered 28 and buys with probability 0.168. Subtracting $\lambda$ says each unit sold today forgoes about 20.71 of later value, so a sale at 16 is worth $16 - 20.71 = -4.71$ net and only higher prices pay. That is how the policy slows sales from 0.595 to 0.104 per arrival.
• The online loop raises $\lambda$ by $\eta(1 - 0.12) = 1.76$ after every sale and lowers it by $\eta \times 0.12 = 0.24$ after every non-sale, so it drifts toward the $\lambda$ at which sales happen at the planned 12% rate. That is the same condition the bisection solved for, and the average $\lambda$ over the second half of the run (21.18) lands next to the bisection value (20.71).
• The online learner earns 3013.6: more than the high-end fixed shadow price (2843.8), less than the low end (3107.4) and the mixture (3090.4), which use the true demand curve. Its $\lambda$ keeps moving between neighbouring prices, which is why it lands between the two ends, and about 2.5 percent below the mixture is the price of learning demand while selling. It does not meet the plan exactly: the stock ran out before period 1000 in all 200 seeds, at period 913.7 on average (earliest 750, latest 952), so it sells a little faster than planned. The comparison depends on this demand curve and grid.

Validating the policy offline, in order

1. Check overlap. IPS needs every price the new policy might choose to have non-zero probability in the logs. Historical prices set by a rule often have none. Run a small randomised price experiment inside a safe band to create usable logs.
2. Estimate revenue and sales rate with IPS and a doubly robust estimator (a demand-model prediction plus an IPS correction of its residual), with clipped weights (weights capped at a maximum so that one row cannot dominate, at the price of some bias), the effective sample size (the number of equally weighted rows that would carry the same information; a small value means a few rows decide the estimate), and confidence intervals.
3. Simulate the inventory path. Off-policy estimates are per-arrival averages; whether stock runs out in week two is a statement about a trajectory, so replay arrivals through a fitted demand model with seeds and compare the stock-out time and the leftover stock.
4. Shadow mode, then a ramp against a holdout that keeps the current pricing.

Guardrails

Price bounds and step limits; a stock-pace alarm that freezes learning and falls back to the safe pricing if sales deviate from plan; a revenue and conversion floor against the holdout, tested sequentially (checked repeatedly as data arrives, with a stopping rule that accounts for the repeated looks); a cap on exploration width; logging of every price with its probability; an audit for unequal outcomes across customer groups; a manual override.

Continuous-action extensions

These are for when a grid of a few prices is too coarse; the grid design above is enough when the useful range is covered by a handful of prices.

• Gaussian-process bandits treat revenue as a smooth function of price and use a kernel (a similarity function that says how strongly evidence at one price informs a nearby price) so that evidence at one price narrows uncertainty at nearby prices; the upper-confidence rule picks the next price. Exact inference scales with the cube of the number of observations, so approximate versions are used at scale.
• Policy gradient (a parametric price policy such as a Gaussian whose mean is a function of context, trained to raise expected reward) fits continuous prices directly. It needs the logging density of each price to be evaluated offline, and the exploration noise is itself a business risk to bound.
• Both need the same inventory multiplier; the dual update does not change.

Common mistakes

Greedy revenue maximisation that ignores stock; evaluating only revenue and not the pace constraint; IPS on logs that never tried the prices being evaluated; letting exploration reach prices outside a trusted band.

Running the code

Use Python 3.12 with numpy 2.1.3. Save each listing to a file and run it with python; every random draw uses a fixed seed.

Direct answer

An epsilon-greedy agent keeps, for each of $k$ arms (the choices it can pull, such as ads or articles), a count of pulls and a running average reward. At each step it flips a coin: with probability $\varepsilon$ (epsilon) it explores by pulling a random arm, otherwise it exploits by pulling the arm with the highest average so far. After seeing the reward it updates only that arm's count and average. Bernoulli rewards mean each pull returns 1 (a success, such as a click) or 0.

Plain-language pseudocode (for readers who do not code)

1. Start every arm with 0 pulls and an average of 0.
2. Draw a random number between 0 and 1. If it is below epsilon, choose any arm at random.
3. Otherwise choose the arm with the highest average (if several tie, choose one of them at random).
4. Pull the arm and observe 1 or 0.
5. Add 1 to that arm's pull count, and move its average toward the reward: new average = old average + (reward - old average) / pulls.
6. Go back to step 2.

Spreadsheet version. Make one row per arm with columns pulls, wins and mean = wins / pulls. Each step, pick a row by the rule above, add 1 to its pulls, add the reward to its wins, and recompute its mean; no other row changes. The program below prints the first six steps this way.

Implementation, with a check against a batch recomputation

The incremental update $\bar x_n = \bar x_{n-1} + (x_n - \bar x_{n-1})/n$ gives the same average as summing all rewards and dividing by $n$, without storing the history. The program verifies that against the logged history. It then runs one 2,000-step simulation and, separately, compares four epsilon values over 100 seeds. True success rates are 0.20, 0.35, 0.50, 0.30 and 0.10, so arm 2 (index 2, the third) is best and a perfect player would earn about $2000 \times 0.5 = 1000$ successes in expectation.

import numpy as np

class EpsilonGreedy:
    def __init__(self, k, epsilon, rng):
        self.k, self.epsilon, self.rng = k, epsilon, rng
        self.counts = np.zeros(k, dtype=int)          # pulls per arm
        self.values = np.zeros(k)                     # running mean reward per arm

    def select(self):
        if self.rng.random() < self.epsilon:
            return int(self.rng.integers(self.k))     # explore: uniformly random arm
        best = np.flatnonzero(self.values == self.values.max())
        return int(self.rng.choice(best))             # exploit: best mean, ties broken at random

    def update(self, arm, reward):
        self.counts[arm] += 1
        self.values[arm] += (reward - self.values[arm]) / self.counts[arm]   # incremental mean

def simulate(true_p, epsilon, steps, seed):
    rng = np.random.default_rng(seed)
    agent = EpsilonGreedy(len(true_p), epsilon, rng)
    history = []
    for _ in range(steps):
        arm = agent.select()
        reward = int(rng.random() < true_p[arm])      # Bernoulli reward
        agent.update(arm, reward)
        history.append((arm, reward))
    return agent, history

true_p = np.array([0.20, 0.35, 0.50, 0.30, 0.10])
agent, history = simulate(true_p, epsilon=0.1, steps=2000, seed=0)
arms = np.array([h[0] for h in history]); rewards = np.array([h[1] for h in history])
print("first 6 steps (arm, reward):", history[:6])
# Spreadsheet view: replay the first 6 steps and show the sheet row that changes
pulls = [0] * 5; wins = [0] * 5
print("step arm reward  pulls_of_arm  wins_of_arm  mean_of_arm")
for step, (a, r) in enumerate(history[:6], start=1):
    pulls[a] += 1; wins[a] += r
    print(f"{step:<4} {a:<3} {r:<6}  {pulls[a]:<12}  {wins[a]:<11}  {wins[a] / pulls[a]:.3f}")
print("pull counts per arm:", agent.counts.tolist(), "sum =", int(agent.counts.sum()))
print("cumulative reward:", int(rewards.sum()))

# Check against a batch recomputation from the history
batch_mean = np.array([rewards[arms == a].mean() if (arms == a).any() else 0.0 for a in range(5)])
print("incremental mean equals batch mean:", bool(np.allclose(agent.values, batch_mean)))
assert (np.bincount(arms, minlength=5) == agent.counts).all()

# Does epsilon matter? Average over 100 seeds of 2000 steps (cumulative reward and share of pulls on the best arm).
print("epsilon  mean cumulative reward  mean share of pulls on best arm")
for eps in [0.0, 0.01, 0.1, 0.3]:
    tot, share = [], []
    for s in range(100):
        a, h = simulate(true_p, eps, 2000, 1000 + s)
        tot.append(sum(r for _, r in h)); share.append(a.counts[2] / 2000)
    print(f"{eps:<7}  {np.mean(tot):8.1f}                {np.mean(share):.3f}")

Output:

first 6 steps (arm, reward): [(2, 1), (1, 0), (2, 0), (2, 0), (2, 0), (1, 1)]
step arm reward  pulls_of_arm  wins_of_arm  mean_of_arm
1    2   1       1             1            1.000
2    1   0       1             0            0.000
3    2   0       2             1            0.500
4    2   0       3             1            0.333
5    2   0       4             1            0.250
6    1   1       2             1            0.500
pull counts per arm: [53, 45, 1814, 42, 46] sum = 2000
cumulative reward: 969
incremental mean equals batch mean: True
epsilon  mean cumulative reward  mean share of pulls on best arm
0.0         690.5                0.310
0.01        836.0                0.581
0.1         925.8                0.832
0.3         862.4                0.731

How to read the output:

• Spreadsheet rows. Step 1 pulls arm 2 (reward 1, so its mean is 1.000). By step 5 arm 2's mean has fallen to 0.250 after four pulls. At step 6 the agent pulls arm 1, whose mean is still 0.000, while arm 2's mean is the highest (0.250). A greedy choice would not do that, so step 6 is an exploration step.
• Pull counts. Arm 2 gets 1,814 of the 2,000 pulls, and the other four arms each get 42 to 53. That is the epsilon-greedy signature: roughly $\varepsilon/k$ of the pulls (here 0.1 / 5 = 2% of 2,000 = 40 on average) go to every arm, and nearly everything else to the current best. The counts sum to 2,000, as the check requires.
• Cumulative reward is 969 for this seed.
• Does epsilon matter? Averaged over 100 seeds, epsilon = 0 earns 690.5 and puts only 31.0% of pulls on the best arm: with no exploration the agent often locks onto whichever arm paid off first. Epsilon = 0.1 earns 925.8 and puts 83.2% there; epsilon = 0.3 earns 862.4 because it wastes more pulls. 0.1 is a common starting value, not a universal default (the class takes epsilon as an argument on purpose); tune it for your horizon and reward rates.

Complexity

• select is $O(k)$ (finding the maximum and the ties); update is $O(1)$. A run of $T$ steps costs $O(Tk)$ time and $O(k)$ memory (the history list in the demo is only for the check). With a heap or a maintained maximum, selection can be made faster, but for small $k$ it is not worth it.

Pitfalls

• Ties. argmax always returns the first index, which biases early pulls toward arm 0. The code breaks ties at random.
• Fixed epsilon has linear regret (regret means expected reward lost compared with always pulling the best arm): exploration never stops. Decay epsilon over time if the arms do not drift.
• Non-stationary rewards. If success rates change, replace the $1/n$ step with a constant step size so old rewards fade.
• Reproducibility. Give the generator a seed, as here, so a bug is repeatable.