Food deliveries, Bayes and Computational Statistics

Introduction

I was grabbing a burger at Shake Shack, Mall of the Emirates in Dubai, when I noticed this notebook on the counter. The staff is using it to track food deliveries, and each service (Carriage, Talabat, UberEats, Deliveroo) has its own column with the order numbers. Let's assume this is the only page for the day, and ask ourselves: given this data, what is the probability that UberEats is the most popular food delivery service?. The ipython notebook is up on Github.

This is good challenge to exercise our Bayesian reasoning. Bayes theorem states that $ P(A|B)P(B) = P(A,B) = P(B|A)P(A) $, where $ P(A|B ) $ is the conditional probability of A given B, and $ P(A, B) $ is the joint probability of A and B co-occuring. All we have to do is replace A with H for Hypothesis and B with D for Data, and we get our formula for attacking the challenge: $ P(H|D) = \frac { P(D|H)P(H) }{ P(D) } $. Here Hypothesis is that UberEats is the most popular service, and Data is the observed counts (2, 2, 16, 14) in the notebook.

Simulation

First, let's do a brute-force simulation. Our model is that there are 4 services, each service has an associated probability, the probabilities sum to 1. When a new orders arrives at Shake Shack, it is according to this model that the service is drawn. A single simulation run---given a model (so, 4 probabilities for each service, summing to 1)---consists of 2 + 2 + 16 + 14 = 34 draws, and the outcome is the distribution between the services, eg. (10, 10, 10, 4). In mathematics terms, this drawing is called a multinomial distribution, numpy has a multinomial function, so we will use that. So our outer loops looks something like:

from random import random
from numpy.random import multinomial

data = {
    'Carriage':  2,
    'Talabat':   2,
    'UberEats':  16,
    'Deliveroo': 12,
}

n_models = 1000

def random_model(k):
    ps = [random() for _ in range(k)]
    s = sum(ps)
    return [p/s for p in ps]

def random_simulation(ps, n):
    return list(multinomial(n, ps))

for i in range(n_models):
    model = random_model(len(data.keys()))
    # call random_simulation(model, sum(data.values())) and count things

Since this is a simulation, we had to introduce n_model, the number of random models we sample (the ensemble size). Now we just have to implement the $ P(H|D) $ formula. This is simple. $ P(D|H) $ is the fraction of random_simulation cases that match the data, when the model matches the hypotehesis (see next sentence). $ P(H) $ is the fraction of models (!) that match the hypothesis of UberEats being the most popular, ie. that probability being the highest. $ P(D) $ is the total number of cases when random_simulation matches the data. First, some helper functions:

def same_as_data(outcome):
    return all([i == j for  (i, j) in zip(outcome, list(data.values()))])

def model_satisfies_hypothesis(ps):
    return max(ps) == ps[2] # UberEats

We have to introduce another parameter, the number of trial runs per model, which we'll call n_simulations. The simulation is now (note that this is actually two loops, the list expression for hits is a loop):

n_models = 1000
n_simulations = 1000*1000

hits_total = 0
hits_hypothesis = 0
num_model_satisfies_hypothesis = 0
for i in range(n_models):
    model = random_model(len(data.keys()))
    hits = sum(
        [same_as_data(random_simulation(model, sum(data.values()))) for _ in range(n_simulations)]
    )
    hits_total += hits
    if model_satisfies_hypothesis(model):
        hits_hypothesis += hits
        num_model_satisfies_hypothesis += 1

In terms of these variables:

p_data_given_hypothesis = hits_hypothesis / (num_model_satisfies_hypothesis * n_simulations)
p_hypothesis = num_model_satisfies_hypothesis / n_models
p_data = hits_total / (n_models * n_simulations)
# p_hypothesis_given_data = p_data_given_hypothesis * p_hypothesis / p_data
# but, for better numerics, the above simplifies to:
p_hypothesis_given_data = hits_hypothesis/hits_total

That's it! We just have to run this simulation, and it will tell us our desired probability. Unfortunately, this is actually a lot of random(), calls. On my Macbook, this takes a long time to run. So I modified it a bit, updated and displayed the estimate hits_hypothesis/hits_total as the simulation was running. This is what it looks like:

Obviously it takes a long time for the simulation to converge. In fact, the p_est shown in the screenshot above is not correct, the simulation has not yet converged. The reason is that the simulation is sparse. In most cases, there are no hits, so we waste a lot of cycles not updating our estimate hits_hypothesis/hits_total.

Monte Carlo integration

The above code is a good start, but it's terribly inefficient. Actually, there's no reason to brute-force the simulation itself, since there is an explicit formula for the multinomial probability of drawing $ k_i $ given $ p_i $ probabilities. The package scipy implements it, it's also called multinomial (like the numpy one), but this is the explicit probability mass function, not a random draw. Actually all we have to do is replace one line, the hits calculation, to compute the fractional probability of observing the data. It changes the meaning of the variables (the hits variables are no longer integers counts, instead they're summed fractions), but the ratios are the same. So the simulation becomes:

from random import random
from scipy.stats import multinomial

data = {
    'Carriage':  2,
    'Talabat':   2,
    'UberEats':  16,
    'Deliveroo': 12,
}

n_models = 100*1000

def random_model(k):
    ps = [random() for _ in range(k)]
    s = sum(ps)
    return [p/s for p in ps]

def same_as_data(outcome):
    return all([i == j for  (i, j) in zip(outcome, list(data.values()))])

def model_satisfies_hypothesis(ps):
    return max(ps) == ps[2] # UberEats

hits_total = 0
hits_hypothesis = 0
num_model_satisfies_hypothesis = 0
for i in range(n_models):
    model = random_model(len(data.keys()))
    hits = multinomial(sum(data.values()), model).pmf(list(data.values())) # <-- this changed
    hits_total += hits
    if model_satisfies_hypothesis(model):
        hits_hypothesis += hits
        num_model_satisfies_hypothesis += 1
p_hypothesis_given_data = hits_hypothesis/hits_total
print(p_hypothesis_given_data)

Note that we now no longer need the variable n_simulations. This is now lighting fast, converges in a few seconds, and quickly tells us that the probability we're looking for is 72%. Given what I saw on the notebook, the probability that UberEats is their most popular delivery service is ~72%. Given this code, it's easy to check another one. What's the probability that Deliveroo is the most popular? Just change model_satisfies_hypothesis, run it again, and it's ~28%. The other two have almost zero chances of actually being the most popular. In this whole discussion I'm assuming a uniform prior, ie. I have no prior bias towards believing one service being more popular than the rest.

Note: this post has a subtle bug in the prior assumption. For more, see this subsequent post. The prior, as implemented above, does not favor one courier over another, but it is not in fact uniform in p-space.

Before moving on, here's a screenshot of this version, but before it converges. It shows the initial oscillations, which is what we were observing with the brute-force simulations (note that it eventually converges to 0.72):

Integrals

In case you remember college calculus, what we're doing above is a Monte Carlo (MC) integration of:

$$ \frac { \int_{\sum p_i = 1 \wedge p_3 = max(p_i)} P(x_i=k_i) dp_i }{ \int_{\sum p_i = 1} P(x_i=k_i) dp_i } $$

where:

$$ P(x_i = k_i) = \frac {N!}{x_1! ... x_s!} p_1^{x_1} ... p_s^{x_s} $$

is the multinomial probability mass function evaluated at

$$ k_i = (2, 2, 16, 14), s = 4, i = 1.. 4, N = \sum k_i = 34.$$

It's called Monte Carlo because instead of computing the integral explicitly, we're taking random samples in p-space and using that to estimate. The $ P() $ function above can be integrated, it's just a polynomial, and the whole integral can be explicitly calculated (I think) with pen and paper (or Mathematica); both numerator and denominator constraints are a hyperplane in 4D p-space, with the numerator being more tricky because of the $ max() $ constraint. I will let the reader do this as homework :)

If you're having trouble understanding this post, check out the free book Think Bayes.