We have a class with \(N=19\) children and one biscuit. We would like to select one child to give the biscuit, but fate should decide which child. Can we use a fair six-sided die to select a winner, with multiple throws, such that all children have an equal probability of winning? And if so, what strategy would be best, in the sense that it minimizes the expected number of throws?

To answer these questions, we start with the analysis of simple cases to learn, hopefully, how to handle the case with \(19\) children. In this process we design a simple and intuitive strategy. As this policy is somehwat hard to analyze mathematically, we'll use simulation to assess its quality. It turns out that this policy is unfair: some children will have a larger probability of winning than others. Then we retract a bit and provide a better algorithm, which is fair and optimal in the sense that is requires the least number of throws on average to select the winner. {{{sidenote{(I discussed this problem with Henk Tijms and he pointed out an answer on Quora on how to handle a class with 20 pupils. Bases on this, I generalized this answer to an algorithm that can deal with an arbitrary number of children.)}}}

In this section we introduce some general concepts that are generally useful. The first is recursion, the second is to use acceptance-rejection policies to decide whether a game should continue or stop.

The first aim is to find a policy that selects in a fair way a winner. But before embarking on how to do this for general classes, we start with considering some simple problems. Perhaps we spot a pattern!

Obviously, when there is \(N=1\) child, there is no need to choose. The next simplest case is when we have \(N=6\) children: number the children from 1 to 6, and throw the die.

When \(N=4\), we can use the following policy. Throw the die for each child. Each child whose outcome is below \(5\) `survives' to the next round. When there are no survivors after the first round (all threw \(5\) or \(6\) ) let them all throw once again. When there are \(1\), \(2\), or \(3\) survivors, use the suitable strategy we developed above. Thus, the idea of this algorithm is like this: use the die to reduce the group such that a set of children remains for which we already have a policy. However, while this is a fool-proof strategy, it requires at least \(4\) throws. We'd better search for more efficient policies.

Note that the first throw splits the group into two subgroups \(\{1,2\}\) and \(\{3, 4\}\) of the same size. Then we have two smaller groups, and for these smaller group we already know how to find a fair solution. Clearly, the pattern of recursion is to try to express the solution in terms of solutions of problems we solved earlier.

This second strategy is known as an acceptance-rejection policy: accept an outcome of an experiment when the outcome lies in a certain subset of the state space, otherwise reject the outcome, and sample again.¹

We can compute the expected number of throws for the second policy with, so-called, first-step analysis. Write \(Y\) for the number of throws we need for the die to hit the set \(\{1,2,3,4\}\). The probability to hit this set is \(p=2/3\), while the probability to hit \(\{5,6\}\) is \(q=1-p\).

In probability terms, we say that \(Y\) has the first success distribution with failure probablity \(q\). This is rv is similar to a geometric rv, but slightly different. The first success rv counts the number of throws up and including the number of throws until a success; the geometric rv counts the number of failures until the success, so one less.² From this,³

So, for \(N=4\), \(\E Y = 3/2\).

Interestingly, the first strategy, based on recursion, always needs exactly \(2\) throws. Instead, the acceptance-rejection method requires less throws on average, but there is no formal bound on \(Y\). We remark in passing an important take away in the design of random algorithms. Sometimes it is required that an algorithm finishes in a guaranteed fixed amount time, rather than that it finishes fast. When time constraints are relevant, it's good to have a deterministic algorithm, even though it can be slower on average.

Observe that with the strategy discussed in 1 for \(N=19\), we always have to throw at least three times. However, the strategy of 1 sometimes completes in just two throws. Might there be a optimal policy that needs somewhere between \(2\) and \(3\) throws on average?

One promising idea to improve on the strategy of using complete rounds of two throws for \(N=19\) is as follows. If there is no winner after the first and second throw, just throw the die a third time. Now check whether the result of the second and third throw is in the list of winning outcomes. To illustrate, if \((X_1, X_2) = (4, 3)\), then, by 1, the outcome \((4,3)\) is not a winning sequence, so that we have to continue with a third throw. After seeing \(X_3=4\), there is winning sequence, since \((X_2, X_3) = (3,4)\), which corresponds to child \((3-1)\cdot 6 + 4 = 16\). If, however, \((X_1, X_2, X_3) = (4, 5, 1)\) there is not yet a winner, because neither \((4,5)\) nor \((5,1)\) correspond to a child. Supposing that \(X_{4}=3\), the result of the last two throws is \((X_3, X_4) = (1, 3)\) so child \(3\) wins.

Notwithstanding that this is an appealing strategy, I was not convinced that this strategy will result in all children having an equal probability to win. The reason is that it reminded me of the, so-called, `monkey-typing problem'. In a later chapter we will explore this problem, but since its analysis requires some new mathematical ideas, we will use here simulation to see how our this new `monkey' policy performs when \(N=19\).

Before we discuss the code below, you should know that in Python (like in C and many other programming languages) arrays start at index 0. This will turn out very practical later in this chapter when we use modulo computions. Hence, in all maths and the code below, children start couting at \(0\), and if we have, for example \(4\) children, the associated set will be \(\{0, 1, 2, 3\}\), and not \(\{1, 2, 3, 4\}\).

The main parts of the next python code work as follows. The tuple last_two corresponds to the outcomes of the last two throws. The rejection set reject contains the tuples that do not correspond to a winner, hence require at least one further throw; hence, reject contains the outcome pairs \((3,1), \ldots, (5,5)\). If last_two is an element of reject, we throw the die again and update last_two to include the last throw. We continue with throwing until last_two is accepted. Once we have a winner, we add one win for this child. Finally, the count object counts how often each outcome occurs.⁴

import collections
import itertools

import numpy as np

rng = np.random.default_rng(1)

reject = set().union((3, i) for i in range(1, 6))
for i, j in itertools.product([4, 5], range(6)):
    reject.add((i, j))

count = collections.Counter()
num_runs = 100_000
for _ in range(num_runs):
    last_two = (rng.integers(6), rng.integers(6))
    while last_two in reject:
        last_two = (last_two[1], rng.integers(6))
    count[last_two] += 1
most = count.most_common()[0][1]
least = count.most_common()[-1][1]
mean = num_runs // 19

print(f"{len(reject)=}, {len(count)=}")
print(f"{least=}, {mean=}, {most=}, {count.total()=}")

len(reject)=17, len(count)=19
least=4464, mean=5263, most=5677, count.total()=100000

We print three checks: the size of the rejection set is \(17\), as it should, the counter counted \(19 = 35 - 17\) different outcomes, and we see \(10^5\) samples in total. The difference between the children that win the most and the least is about 900, i.e., about \(1\%\) for a sample size of \(10^{5}\) throws. Is this large? To provide some context, we can compare this with a simulation of uniform numbers on \(\{0, 1, \ldots, 18\}\).

import numpy as np

rng = np.random.default_rng(1)
N = 100_000
X = rng.integers(19, size=N)
bins = np.bincount(X)
print(f"{min(bins)=}, {max(bins)=}")

min(bins)=5144, max(bins)=5363

This variation between the numbers that occur the most and least often is quite a bit smaller than the variation under the monkey strategy. We can do further statistical tests, but for the moment I simply conclude that, under the monkey-typing strategy, the children don't have uniformly distributed winning chances.

Still I find it interesting to analyze how many throws this monkey strategy needs on average tof inish, in particular how it compares to the acceptance-rejection policy designed in 1. If it slower, then the policy to use (multiple) rounds of three throws times is better in all respects, but if the monkey policy is faster, we might make a trade off between fairness and speed. Note that this is a common theme in optimization: nearly always there are multiple objectives that we like to minimize or maximize, but we cannot perform optimally on all dimenstions at the same time. If so, we need to make a choice which objective we prefer over another.⁵

Returning to the analysis of the average number of throws it takes the monkey to select a winner, consider the code below. The next exercise asks you to explain some parts of it.

T = 0
for _ in range(N):
    last_two = (rng.integers(6), rng.integers(6))
    T += 2
    while last_two in reject:
        last_two = (last_two[1], rng.integers(6))
        T += 1

print(f"{T/N=}")

T/N=2.94335

As T counts the total number of throws for \(N\) games, T/N is the average number of throws per game. Apparently, this lies just a bit below \(3\). In summary, although unfair, the monkey policy needs less throws on average than the best acceptance-rejection policy so far, as this already needs three throws for the first round.

In the above we learned that when the class contains \(N=19\) children, the monkey strategy, which prescribes to continue throwing the die until a winning pattern emerges, does not give a fair outcome. Hence, for the moment, by 1 the best fair algorithm for \(N=19\) seems to throw three times: if that gives a winner, stop; otherwise throw three times again, and continue with rounds until there is a winner.

However, there is a strategy that is fair and minimizes the expected number of throws at the same time. This strategy we discuss now. Although a bit harder to understand, it's an elegant interplay between probability and number theory.

In the sequel we write \(Z\sim\Unif{\{0, 1, \ldots, n\}}\) to say that a rv \(Z\) is uniformly distributed on the set \(\{0, 1, \ldots, n\}\). We also write \(Z|A\) to denote a rv \(Z\) given that the event \(A\) materialized.

Let us first estimate the number of throws needed for a class with \(N=20\) children; this is somewhat simpler than the case with \(N=19\). Consider the following code, where \(X\) is an iid uniform rv on \(\{0,\ldots,5\}\) corresponding to the throw of the die. The line t % 20 computes the remainder of t after dividing by \(20\).

import numpy as np

rng = np.random.default_rng(1)

t = 0
while t < 16:
    X = rng.integers(6)
    t = 6 * t + X
t = t % 20

The above two exercises explain the threshold \(16\) of the algorithm for \(N=20\). Now that we understand this, we can try to generalize the algorithm to other values for \(N\), such \(N=19\). Consider the following update rule for \(i=1, 2, \ldots\), \(t_0=0\), and \(X_1\sim\Unif{\{0, \ldots, 5\}}\),

Notice that for \(N=19\) the threshold is not the same for each throw, i.e., 17 for the second and 7 for the third throw.

Clearly, to deal with cases such as \(N=19\) we need to modify the above algorithm. Here is a generalized version.⁶

def select_uniformly_with_a_die(N):
    n = r = 1
    t = 0
    while t < r:
        X = rng.integers(6)
        t = 6 * t + X
        n = 6 * r
        r = n % N
    return t % N

Let us test this generalized algorithm numerically. Here is the code.

count = collections.Counter()
for _ in range(num_runs):
    s = select_uniformly_with_a_die(19)
    count[s] += 1

most = count.most_common()[0][1]
least = count.most_common()[-1][1]
mean = num_runs // 19

print(f"{len(count)=}, {count.total()=}")
print(f"{least=}, {mean=}, {most=}")

len(count)=19, count.total()=100000
least=5144, mean=5263, most=5434

If we compare this output to our earlier results, the difference between the least and most winning child is reasonable. In fact, comparing the occurrences, we see that our algorithm performs very well!

The last step is to estimate the expected number of throws \(\E Y\) required for the case \(N=19\). First, however, we consider \(N=20\) because in this case we can obtain a simple closed-form expression for \(\E Y\).

To compute \(\E Y\) for general \(N\) we can slightly modify the above algorithm.

N = 20
n = r = 1
alpha, eps = 1, 1e-3
EY = 0
while alpha >= eps:
    EY += alpha
    n = 6 * r
    r = n % N
    alpha *= r / n
    print(f"{n=:3d}, {r=:3d}, {r/n=:2.4f}, {alpha=:2.4f}, {EY=:2.4f}")

n=  6, r=  6, r/n=1.0000, alpha=1.0000, EY=1.0000
n= 36, r= 16, r/n=0.4444, alpha=0.4444, EY=2.0000
n= 96, r= 16, r/n=0.1667, alpha=0.0741, EY=2.4444
n= 96, r= 16, r/n=0.1667, alpha=0.0123, EY=2.5185
n= 96, r= 16, r/n=0.1667, alpha=0.0021, EY=2.5309
n= 96, r= 16, r/n=0.1667, alpha=0.0003, EY=2.5329

This seems to converges to \(2.533\). The last task is to run the code for \(N=19\). This gives \(\E Y=2.509\). Comparing this to the monkey strategy of 1 we see that this algorithm does better on all respects: it is fair and requires less throws on average.

And now two final questions.

By studying the above toy problem we covered many nice and useful ideas. We developed acceptance-rejection strategies to select a winner. Then we tried to speed up this algorithm to the monkey algorithm, but this algorithm turned out to be unfair. However, by improving the idea, we were able to find a fair and optimal algorithm. This final algorithm is not quite trivial, but effective and suprisingly short.

In the next section we'll look at the monkey-typing strategy. This will provide us with a number of useful tools, in particular first-step analysis.

In a number of my org mode files I include exercise with hints and solutions. In the pdf files, which I make via \(\LaTeX\), I use the answers package print the exercises in the main text and the hints and solutions at the end of the file. However, when exporting the file to html via my website I would like to see the exercises, hints and solutions printed like so:

I use emacs and org mode to write files, and nikola to make my homepage. After some searching on the web, I found the org-special-block-extras package which proved promissing for my goal. The related github page shows how to install it. Here are the steps I followed.

First I installed the org mode plugin for nikola. The orgmode.py file shows that the next command can used to convert an org to an html file that nikola can read.

/usr/bin/emacs --batch -l init.el --eval '(nikola-html-export "input.org" "output.html")'

The org mode plugin also provides an init.el file that I used as a starting point.

Since I need functionality from org-special-block-extras, I have to update the init.el anyway, so here it is.

Let's test this on this file. As this file sits in the posts we need to include the path to init.el.

/usr/bin/emacs --batch -l ../plugins/orgmode/init.el --eval '(nikola-html-export "making-an-exercise-environment.org" "../output/output.html")'

I have a house for sale and I am prepared to organize at most \(n\geq 1\) house viewings. Assume that offers are idd rvs \(\sim \Unif{[0,1]}\), and that I have to accept an offer on the spot or reject it. If rejected, the offer is retracted, hence I cannot accept it at a later date. Under these rules, what is the best price I can expect?

Once my house is sold, I want to buy a new house. As this is a multi-objective optimization problem, I don't want to use just the price of the house to decide to buy it or not. In this case I just assume I can rank the houses from best to worse. Thus, I cannot use the same policy as when selling a house, but need some other policy to decide.

Both these problems can be solved as an optimal stopping problem. So, let's see how to do this, partly with Snell's envelope. I include some python code to get numerical output.

When selling a house, at most \(n\) buyers will present their offers \(\{X_i\}_{i=1}^n\) , iid and uniform on \([0,1]\), one after the other. Rejected offers are lost. What is the best price? I used this site, where the same problem is discussed in terms of the Secretary Problem.

Contrary to when selling a house, when buying a house I don't have an absolute benchmark. In the house selling case, all offers are capped by \(1\), and when the very first bid is already very close to \(1\), I know that that the probability of getting a better offer at a later date will be small. However, when buying a house, I assume I can just rank the houses, but I don't have benchmark. Therefore I will use the first \(r\) say visits out of maximally \(n\) to sample the quality of the houses.

To understand how to apply Snell's envelope to the house selection problem, I used these two resources besides `Promenade aléatoire' by M. Benaïm and N. El Karoui.

• Bodineau, Modélisation de phénomènes aléatoires
• Ferguson lecture notes. These notes provide a formal statement of the house selection problem and some interesting discussions on further work.

These accounts are clear, but leave quite a bit of work to the reader. The idea here is to also do this work.

We assume that \(n > 2\), because when \(n=1,2\) the problem is not interesting.

A Vickrey auction is a sealed-bid auction where bidders submit bids without knowing the bids of other people. The winner is the participant with the highest bid, but interestingly, the price paid is the second-highest bid price rather than the highest. One way to understand this is to think about an auction in which bidders accept a price announced by the auctioneer by raising their hand. If the second highest bidder is prepared to pay \(A\), then s/he will stop when the price is \(A+1\), and this is the price at which the highest bidder gets the item. Clearly, the highest bidder does not reveal its own maximum price to win the auction; in fact, all bidders except the highest reveal their bid.¹

If we have large debts, and an expensive car, should we organize a Vickrey auction to sell our car and pay our debts? How will the expected value and the variance of the second bid depend on the number of bidders? Perhaps we need a minimum number of bidders to reduce the probability of getting a low value for our car. To obtain insight into these questions we will make a histogram of simulated prices at which the car is sold in a Vickrey auction.

As a simple model suppose we have \(10\) bidders whose bids are idd rvs \(\{X_i\}_{i=1}^{10}\) according to some distribution \(F\), say. To estimate the value of the second largest bid by simulation, we make a matrix \(X\) of size \(N\times 10\) with entries \(X_{i j} \sim F\); like this, each row corresponds to one round of an auction, and column \(j\) to the bids of the \(j\)th participant. Finding the hightest bid per row is easy with the code np.max(X, axis=1), but for the Vickrey auction we don't need the largest, but the second largest element of each row. One way to find this is to sort each row, and then pick the values of \(9\)th column of \(X\), i.e., X[:,8]. However, this is algorithmically more expensive than necessary: we don't need to sort all values of the row, we only need to the know the, so-called, \(k\)th median, with \(k=9\).

To build a better understanding of how to simulate efficiently the Vickrey auction, we first need to discuss some elementary, but highly useful, theory to express the complexity of algorithms. Then we apply this theory to an algorithm to solve the Traveling Salesman Problem (TSP) by means of recursion. As we will see, the implementation of the TSP is really elegant. Then we consider the complexity of two sorting algorithms. The second of these, quicksort, provides us with a very nice and efficient recursive algorith to find the \(k\)th median in an array of numbers. In the last section we use this algorithm to find the second highest bid in a Vickrey auction and plot of the density of these bids to obtain insight into the price we can expect for our car.

In simulation we often repeat the same algorithm with different so that we obtain statistical insight by means of many sample paths. Therefore the difference in time required by a dumb and a smart algorithm to complete one simulation run can be huge. One way to compare the speed of algorithms is by means of the concept of complexity, for which we use the order symbol \(O(\cdot)\). Here are two examples to illustrate the idea. First, when adding two \(n\)-digit numbers, we add the \(m\)th digit in the first number with the \(m\)th digit in the second number. We do the same operation for all digits, and since there are \(n\)-digits in total, the algorithm requires a number of additions that is \(O(n)\). Second, to multiply two \(n\)-digit numbers by hand, we multiply each digit of the first number by each digit of the second. Clearly this requires \(O(n^2)\) multiplications. Thus, with the \(O(\cdot)\) symbol, we express how the dominating term scales as a function of the input of the data, but we abstract away from the time each specific operation takes. In other words, whether an algorithm has \(4 n^3\) or \(5n^3\) complexity, in all such cases we say it is \(O(n^3)\). It is evident that the larger the power of the \(n\), the more execution time an algorithm need in general to complete for large input.

To practice a bit more with the \(O\) symbol, given a list of \(n\) elements, the following code finds the largest element:

X = [5, 8, 3]
M = X[0]
for x in X[1:]:
    M = max(M, x)

A somewhat more interesting example is the complexity of the traveling salesman problem (TSP). This problem consists of finding the shortest possible route that visits, exactly once, each city of a given set of \(n\) cities, starting from a depot labeled \(0\). The next code compute recursively the length of the shortest path. We start at the depot and provide a set of cities. When in the current city, we solve the tsp that remains after going to each city in cities. (Take your time to understand the recursion.) We add typing to clarity what sort of argument TSP is supposed to accept. By the way, it might seem pointless to overload integers to city_id, but the meaning of the word city_id is much easier to understand.

city_id = int
length = float


def dist(i: city_id, j: city_id) -> length:
    # An arbitrary distance function between cities i and j
    return abs(i - j)


def tsp(current: city_id, cities: set[city_id]) -> length:
    if not cities:  # There are no remaining cities, so return home
        return dist(current, 0)
    L = []
    for city in cities:
        d = dist(current, city) + tsp(city, cities.difference([city]))
        L.append(d)
    return min(L)

depot = 0
cities = set([1, 2, 3, 4])
print(tsp(0, cities))

8

As an aside, the above code can be written as one list comprehension, like so:

def tsp(s: city_id, C: set[city_id]) -> length:
    if not C:
        return dist(s, 0)
    return min(dist(s, c) + tsp(c, C.difference([c])) for c in C)

The line if not cities: requires attention: it is the condition to tell the recursion to stop. More generally, when applying recursion, it is essential to ensure that eventually some stopping criterion is satisfied, for otherwise the algorithm might never stop running. Hence, when desiging recursive algorithms, always check that the algorithm contains a stopping condition.

This closes our discussion of complexity theory. Much more can be said about it, but for our purposes this suffices. We now have the means to compare the efficiency of algorithms.

After the above excursions to identify an efficient algorithm that finds the second highest number \(Y_{2}\), say in a list, we return to our initial questions about selling our car by means of a Vickrey-auction. The aim is to obtain a histogram of the density of \(Y_{2}\)}. Before actually doing the simulation, let us consider the problem from a somewhat higher perspective. Actually, we are interested in the density \(f_{2}\) of the second highest element of the order statistic of the \(\{X_{i}\}\).

In the next code we compute \(f_2\), carry out a simulation and plot the results. We'll use the normal and uniform distributions of scipy.

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm, uniform

from latex_figures import fig_in_latex_format

fig, (ax1, ax2) = plt.subplots(1, 2, sharey="all", figsize=(6, 3))

The parameters are like this. The grid is used for plotting.

bidders = 10
k = bidders - 2
num_sim = 50
gen = np.random.default_rng(3)
grid = np.linspace(0, 1.5, num=50)

We start with assuming that \(X_i\sim\Unif{[0,1]}\). We need the cdf, sf and pdf for the computation of \(f_2\) by means of the formula of the previous exercise.

X = uniform(0, 1)
f, F, G = X.pdf(grid), X.cdf(grid), X.sf(grid)
f2 = 10 * 9 * f * F**8 * G

To simulate, we generate a matrix with bids such that each of its \(N\) rows corresponds to an auction of \(n=10\) bidders, and each element is a bid uniform on \([0,1]\). To get the (simulated) value for \(Y_{2}\), we use the numpy function partition, which uses a fast algorithm to compute the \(k\)th median. The \(k\)-median is at the \(k\)th column of each row, hence the keyword axis=1.

bids = X.rvs((num_sim, bidders), random_state=gen)
sim2 = np.partition(bids, k, axis=1)[:, k]

It remains to make plots.

ax1.set_title("Uniform bids")
ax1.plot(grid, f2, color="black", label="$f_2$")
ax1.hist(sim2, bins=grid, density=True, color="gray", label="sim")
ax1.legend(loc="upper left")

We run virtually the same code for \(X_i\sim \Norm{1/2, 1/\sqrt{12}}\), and plot in Fig. 1 the density \(f_2\) and the histogram obtained by simulation. The right panel of the figure shows the histogram for a simulation with \(X_i\sim\Norm{\mu=0.5, \sigma=1/\sqrt{12}}\). Clearly, the model of the bid sizes has a major effect on the price we can get for our expensive car.

In this section we studied a bit of complexity theory. This allowed us to understand why the quicksort algorithm, based on recursion, is a much better algorithm than bubblesort. Recursion helped us also to develop a nice and efficient method to compute the $k$-median, and to write the TSP in just a few lines. We applied all this to estimating the density of the price of the second bid in an Vickrey auction, because this is the value we obtain for good sold under this auction mechanism. Finally, we realized that the second bid is actually the second value of the order statistic of a number of iid bids.

Actually, I did not expect that the density would have such a big tail to the left, as is shown in Fig. 1. Perhaps this explains why, instead of Vickrey auctions, sellers might prefer Dutch auctions in which the selling price starts very high and decreases until the first person accepts the price. However, let's not draw simplistic conclusions; we should turn to a serious book on auction theory instead. In our work we are occupied with showing the many uses of recursion, some maths, and coding.

Consider the following puzzle: Chose numbers 1 to 5 such that in the diagram below each number is contained exactly once in each row and each column. The numbers in the fields should respect any inequality between two fields. We number the fields from the left top, e.g., field \((1,5)\) contains the number 4.

We will use pulp to solve this puzzle. The sudoku problem discussed in the pulp documentation provides some inspiration on how to tackle our problem.

from pulp import lpSum, LpVariable, LpMinimize
from pulp import LpProblem, LpStatus, value, LpInteger

prob = LpProblem("Number-Puzzle-Problem", LpMinimize)

As decision variables we use \(x_{kij}\) such that \(x_{kij} = 1\) if value \(k\) appears in the field with row \(i\) and column \(j\). The values, rows and columns are all elements of the set \(\{1, 2, 3, 4, 5\}\).

Vals = Rows = Cols =  range(1, 6)
x = LpVariable.dicts("x", (Vals, Rows, Cols), 0, 1, LpInteger)

The simplest constraints deal with the numbers that are given in the diagram.

prob += x[4][1][5] == 1
prob += x[4][4][1] == 1
prob += x[2][4][5] == 1

There should be precisely one value in each field.

for r in Rows:
    for c in Cols:
        prob += lpSum([x[k][r][c] for k in Vals]) == 1

There can be only one value per row.

for k in Vals:
    for r in Rows:
        prob += lpSum([x[k][r][c] for c in Cols]) == 1

And also only one value per column.

for k in Vals:
    for c in Cols:
        prob += lpSum([x[k][r][c] for r in Rows]) == 1

The inequality constraints of the board require more work.

Let's consider an example. In the problem diagram above we see that the value of field \((2,1)\) must be larger than the value in field \((1,1)\). Thus, if field \((2,1)\) has value 3, then field \((1,1)\) is not allowed to have a value of 3, 4, or 5. This means that, in terms of the decision variables, if \(x_{3,21}=1\) then \(\sum_{w=3}^5 x_{w,11} = 0\). More generally, we want that

We can implement this implication with the big \(M\) trick. Choose some big \(M\), then set

so that if \(x_{k, 21} = 1\) then \(\sum_{w=k}^5 x_{w, 11}\) must be zero, while if \(x_{k, 21} = 0\), then \(\sum_{w=k}^5 x_{w, 11}\) is unconstrained. In code this becomes:

M = 100  # This is large enough

for k in Vals:
    prob += lpSum(x[w][1][1] for w in R[k:]) <= M * (1 - x[k][2][1])

In a similar way we implement the other inequalities.

for k in Vals:
    prob += lpSum(x[w][1][5] for w in R[k:]) <= M * (1 - x[k][1][4])
    prob += lpSum(x[w][2][5] for w in R[k:]) <= M * (1 - x[k][1][5])
    prob += lpSum(x[w][4][4] for w in R[k:]) <= M * (1 - x[k][3][4])
    prob += lpSum(x[w][4][1] for w in R[k:]) <= M * (1 - x[k][5][1])
    prob += lpSum(x[w][4][3] for w in R[k:]) <= M * (1 - x[k][4][4])
    prob += lpSum(x[w][4][5] for w in R[k:]) <= M * (1 - x[k][4][4])
    prob += lpSum(x[w][4][1] for w in R[k:]) <= M * (1 - x[k][5][1])
    prob += lpSum(x[w][5][5] for w in R[k:]) <= M * (1 - x[k][4][5])

A much simpler way (which occurred to me quite a bit later) is to define supplementary variables \(v\) that correspond to the value of the cell, i.e.,

Then we must have, for instance, for cells \((1,4)\) and \((1,5)\) that \(v_{14}>v_{15}\). In code this becomes

v = {}
for r in Rows:
    for c in Cols:
        v[(r, c)] = lpSum([k * x[k][r][c] for k in Vals])

prob += v[(1, 5)] <= v[(1, 4)]
prob += v[(1, 1)] <= v[(2, 1)]
prob += v[(2, 5)] <= v[(1, 5)]
prob += v[(4, 4)] <= v[(3, 4)]
prob += v[(4, 1)] <= v[(5, 1)]
prob += v[(4, 3)] <= v[(4, 4)]
prob += v[(4, 5)] <= v[(4, 4)]
prob += v[(4, 1)] <= v[(5, 1)]
prob += v[(5, 5)] <= v[(4, 5)]

Note that since all cells must have different values, all strict inequalities in the constraints can be implemented as a \(\leq\) constraints.

The objective is trivial, because there is nothing to optimize. We only use pulp to find a feasible solution.

prob += 0, "Arbitrary Objective Function"

Solving the problem comes down to calling the solve() method of pulp.

prob.solve()

Have we found a solution?

print("Status: ", LpStatus[prob.status])

Status:  Optimal

Yes! Let's print it.

for r in Rows:
    res = ""
    for c in Cols:
        for k in Vals:
            if value(x[k][r][c]) == 1:
                res += str(k) + " "
    print(res)

1 2 3 5 4
2 4 5 1 3
3 1 2 4 5
4 5 1 3 2
5 3 4 2 1

For some application I needed a random shuffle \(r\) of a list of numbers \(a\). The shuffle, however, should be weighted, such that, given a list of weights \(w\), the element \(a_i\) with weight \(w_i\) should be the first element of the random shuffle \(r\) with probability \(w_i/\sum_k w_k\). I searched for a python implementation of such an algorithm, but I couldn't find one. Also, numpy doesnt' provide it, and numpy.random.shuffle only offers shuffling with uniform weights. Thus, this code.

Weighted random shuffling is the same as weighted random sampling from a list \(a\) without replacement. In steps:

1. Choose with probability \(w_{i}/\sum_{i} w_{i}\) element \(a_{i}\) from \(a\).
2. Add this element to a list \(r\).
3. Remove element \(a_{i}\) from \(a\) and \(w_{i}\) from \(w\).
4. Stop if \(a\) is empty, else,
5. Return to step 1 with the updated lists \(a\) \(w\).

However, I found the code below for weighted random sampling on this page.

from random import random

def weighted_choice(weights):
    rnd = random() * sum(weights)
    for i, w in enumerate(weights):
        rnd -= w
        if rnd < 0:
            return i

Thus, weighted_choice() gives an index to the chosen element of an array. With this index I can fill a random shuffle \(r\) with elements of \(a\).

As removing elements from lists is a bit expensive, I chose to set the weight \(w_{i}\) of the chosen element to zero. This ensures that the element \(i\) cannot be selected the next time.

def weighted_shuffle(a,w):
    w = list(w) # make a copy of w
    if len(a) != len(w):
        print("weighted_shuffle: Lenghts of lists don't match.")
        return

    r = [0]*len(a) # contains the random shuffle
    for i in range(len(a)):
        j = weighted_choice(w)
        r[i]=a[j]
        w[j] = 0
    return r

Once I had the above code, I went to Gerlach van der Heide. He managed to speed it up by several orders of magnitude.

from random import random
from bisect import bisect_right
import numpy as np

def weighted_shuffle(a,w):
    r = np.empty_like(a)
    cumWeights = np.cumsum(w)
    for i in range(len(a)):
         rnd = random() * cumWeights[-1]
         j = bisect_right(cumWeights,rnd)
         #j = np.searchsorted(cumWeights, rnd, side='right')
         r[i]=a[j]
         cumWeights[j:] -= w[j]
    return r

a = np.arange(1,1000)
w = np.arange(1,1000)
r = weighted_shuffle(a,w)

Here is the result

print(r[:2])

Some testing shows, for a reason unknown to me, that np.searchsorted is slower than bisect_right.

While searching the web I came across the one-liner below. It is elegant, but wrong. The weights according to which the elements of \(a\) are chosen are not in proportion to \(w\). Some basic probability theory shows why this is the case. It comes down to proving that \[P\{w_1 U_1 \leq w_2 U_2\} \neq \frac{w_1}{w_1+w_1},\] where \(w_i\) are the weights and \(U_i\) uniform random deviates on \([0,1]\).

def wrong(a,w):
    aa = range(len(a))
    aa.sort(key=lambda i: w[i]*random()*beta[i])
    return [a[i] for i in aa]

We like to compute the number of possibilities \(\P{n, N}\) for \(x = (x_1, \ldots x_n)\) such that \(x_i \in \{0,1,\ldots, N\}\) and \(\sum_i x_i \leq N\). It is easy to see that \(\P{n, N}\) satisfies the recursion: \[\P{n, N} = \sum_{i=0}^N \P{n-1, N-i},\] with boundary conditions \(\P{1, N} = N+1\) for all \(N\). Note that by using the summation above, this condition can be replaced by \(\P{0, N} = 1\) for all \(N\).

Computing \(\P{n,N}\) is easy when we use memoization. In fact, we can code the compuation in nearly one line!

from functools import lru_cache

@lru_cache(maxsize=128)
def P(n, N):
   return (n==0) or sum( P(n-1, N-i) for i in range(N+1) )

n=5
N=80
print(P(n = 5, N = 80))

32801517

We throw multiple times with a coin that lands heads with probability \(p\). What is the probability \(\P{n,k}\) to see at least \(k\) heads in row when you throw a coin \(n\) times?

A bit of thinking shows that \(\P{n,k}\) must satisfy the recursion \[\P{n,k} = p^k + \sum_{i=1}^k p^{i-1} q\, \P{n-i,k},\] because it is possible to throw \(k\) times heads from the first throw, but otherwise you throw \(i\), \(i < k\), times a heads, then a tails, after which you have to start all over again.

Reasoning similarly, the expected number times \(\E{n,k}\) to see at least \(k\) heads in row when you throw a coin \(n\) times must satisfy the recursion \[\E{n,k} = p^k(1+\E{n-k,k}) + \sum_{i=1}^k p^{i-1} q\, \E{n-i,k}.\]

p = 0.5
q = 1. - p

@lru_cache(maxsize=128)
def P(n, k):
    """
    probability to see at least k heads in row when a coin is thrown n times
    """
    if n < k:
        return 0
    else:
        return sum(P(n-i,k) * p**(i-1) * q for i in range(1,k+1)) + p**k

@lru_cache(maxsize=128)
def E(n, k):
    """
    expected number of times to see at least k heads in row when a coin is thrown n times
    """
    if n < k:
        return 0
    else:
        tot = sum(E(n-i,k) * p**(i-1) *q for i in range(1,k+1))
        tot += p**k * (1 + E(n-k,k))
        return tot



k = 2

for n in range(k,10):
    print(n, P(n,k), E(n,k))

Sometimes I want to use a python script to produce a plot and include this plot in \(\LaTeX\). My usual approach was like this. First I wrote a python script that, after some computations, exported the plot to a pdf file (or some other format). Then I imported the pdf file within a figure environment in a \(\LaTeX\) file. Thus, I had three files, of which two were kind of superfluous, i.e., the python script and the pdf file with the figure. And, after some time, I typically forgot which python script I used to make which pdf file, so I had to include comments in my \(\LaTeX\) file to explain where to find what, and what did what.

And then there came a day that I did not like this anymore.

I want one file .tex file that contains all: obviously the story itself, but also the python and the figure, provided of course it takes not a lot of time to make the data for the figure. After some searching and testing, here is an MWE of how I am doing things now with pythontex, matplotlib, and tikzplotlib.

I use arara to compile the LaTeX file, but this just handy, not necessary for the rest to work.

I want to write posts and a homepage in org mode. (I dislike markdown, and I detest typing math in rst.) After watching one of the blogs on Cest la Z I considered nikola as a good option to set things up. It turned out to be really easy; I just followed the steps as described here. I had to install some python packages, but with pip this was a no brainer.

There was one caveat with including source code blocks in org files. After reading the error messages, I noticed that I had to install the emacs htlmize package. The init.el and myinit.org files on my emacs repo show how to set this up.

And now all works!