A problem in probability

20221214-secretary-game.tex

\documentclass[12pt,a4paper]{article}  % -*- coding: utf-8 -*-
\usepackage[a4paper,margin=1.5cm]{geometry}
\usepackage[english]{babel}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{times}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsthm}
%
\usepackage{mathrsfs}
%\usepackage{bm}
%\usepackage{stmaryrd}
\usepackage{wasysym}
\usepackage{url}
\usepackage{graphicx}
\usepackage[usenames,dvipsnames]{xcolor}
%\usepackage{tikz}
%\usetikzlibrary{matrix,arrows,decorations.markings}
%\usepackage{hyperref}
%
%
%
\mathchardef\emdash="07C\relax
\mathchardef\hyphen="02D\relax
\DeclareUnicodeCharacter{00A0}{~}
%
\newcommand{\erf}{\operatorname{erf}}
%
%
%
\begin{document}
\pretolerance=8000
\tolerance=50000

Let $F$ be a fixed probability distribution on $\mathbb{R}$ (we
interpret $F(x)$ as $\mathbb{P}(X\leq x)$).  For $1
\leq k \leq n$, we consider the following one-player
game $\mathbf{S}(n,k)$: the player will be presented with $n$
i.i.d. values $X_0,\ldots,X_{n-1}$ drawn from the distribution $F$; at
each turn $i$, they must decide (knowing $X_i$ and all previous values
and decisions) to either “keep” or “discard” the value $X_i$, in such
a way that, at the end exactly $k$ of the $n$ values have been kept
(and $n-k$ discarded): the final score is the sum of the $X_i$ which
have been kept, and the goal is to maximize this score.

Clearly, in $\mathbf{S}(n,0)$ there is no choice but to discard
whatever values are presented, and in $\mathbf{S}(n,n)$ there is no
choice but to keep.  In $\mathbf{S}(n,k)$ for $0<k<n$, on the other
hand, both options exist, and if the player chooses to discard the
first value $X_0$ they are left playing the game $\mathbf{S}(n-1,k)$
whereas if they choose to keep they are left playing the game
$\mathbf{S}(n-1,k-1)$ with an additional score of $X_0$.

An optimal strategy for $\mathbf{S}(n,k)$ is therefore easily devised
by induction on $n$: if $k=0$ then discard, if $k=n$ then keep, and if
$0<k<n$ then keep $X_0$ when it is greater than a threshold value
$\mathcal{E}(n-1,k) - \mathcal{E}(n-1,k-1)$ where $\mathcal{E}(n,k)$
is the expected value of $\mathbf{S}(n,k)$ for the strategy we are
defining.  (In other words, keep $X_0$ when the expected value $X_0 +
\mathcal{E}(n-1,k-1)$ is better than that $\mathcal{E}(n-1,k)$ we get
by discarding it.)  It is fairly obvious (again, by induction on $n$)
that this strategy is indeed optimal.

Now let us call $G(t)$ the expected value of $\max(X,t)$ when $X$ is
drawn from the probability distribution $F$, in other words
\[
\begin{aligned}
G(t) &= \mathbb{P}(X\leq t)\,t
\;+\; \mathbb{P}(X>t)\,\mathbb{E}(X\,|\,X>t)\\
&= t\,F(t) + \int_{x=t}^{+\infty} x\,dF(x)\\
\end{aligned}
\]
Then $\mathcal{E}(n,k) = \mathbb{P}(X\leq t)\,\mathcal{E}(n-1,k) \;+\;
\mathbb{P}(X>t)\,\Big(\mathbb{E}(X\,|\,X>t) +
\mathcal{E}(n-1,k-1)\Big)$, where $t = \mathcal{E}(n-1,k) -
\mathcal{E}(n-1,k-1)$, is rewritten as $\mathcal{E}(n-1,k-1) + G(t)$.

So we have obtained recurrence relations for $\mathcal{E}(n,k)$,
namely:
\[
\begin{aligned}
\mathcal{E}(n,0) &= 0\\
\mathcal{E}(n,n) &= \mathcal{E}(n-1,n-1) \;+\; G(-\infty)\\
\mathcal{E}(n,k) &= \mathcal{E}(n-1,k-1) \;+\; G\Big(\mathcal{E}(n-1,k)-\mathcal{E}(n-1,k-1)\Big)\\
&\hphantom{=}\text{\quad for $0<k<n$}\\
\end{aligned}
\]
where $G(-\infty) := \mathbb{E}(X)$.  For example:
\[
\begin{array}{c}
\mathcal{E}(0,0) = 0\\
\mathcal{E}(1,0) = 0\;\;\;,\;\;\; \mathcal{E}(1,1) = G(-\infty)\\
\mathcal{E}(2,0) = 0\;\;\;,\;\;\; \mathcal{E}(2,1) = G(G(-\infty))
\;\;\;,\;\;\; \mathcal{E}(2,2) = 2G(-\infty)\\
\mathcal{E}(3,0) = 0\;\;\;,\;\;\; \mathcal{E}(3,1) = G(G(G(-\infty)))
\;\;\;,\;\;\; \mathcal{E}(3,2) = \text{below}
\;\;\;,\;\;\; \mathcal{E}(3,3) = 3G(-\infty)\\
\mathcal{E}(3,2) = G(G(-\infty))+G(2G(-\infty)-G(G(-\infty)))\\
\end{array}
\]

If $F(x) = \frac{1}{2}(1+\erf(\frac{x}{\sqrt{2}}))$ defines a standard
Gaussian distribution ($\frac{dF}{dx} =
\frac{1}{\sqrt{2\pi}}\exp(-\frac{x^2}{2})$), then the expected value
of $\max(X,t)$ is:
\[
G(t) = \frac{t}{2}(1+\erf(\frac{t}{\sqrt{2}})) + \frac{1}{\sqrt{2\pi}}\exp(-\frac{t^2}{2})
\]
Numerically:
\[
\begin{array}{c}
\mathcal{E}(0,0) = 0\\
\mathcal{E}(1,0) = 0\;\;\;,\;\;\; \mathcal{E}(1,1) = 0\\
\mathcal{E}(2,0) = 0\;\;\;,\;\;\; \mathcal{E}(2,1) \approx 0.39894
\;\;\;,\;\;\; \mathcal{E}(2,2) = 0\\
\mathcal{E}(3,0) = 0\;\;\;,\;\;\; \mathcal{E}(3,1) \approx 0.62975
\;\;\;,\;\;\; \mathcal{E}(3,2) \approx 0.62975
\;\;\;,\;\;\; \mathcal{E}(3,3) = 0\\
\mathcal{E}(4,0) = 0\;\;\;,\;\;\; \mathcal{E}(4,1) \approx 0.79041
\;\;\;,\;\;\; \mathcal{E}(4,2) \approx 1.02869
\;\;\;,\;\;\; \mathcal{E}(4,3) \approx 0.79041
\;\;\;,\;\;\; \mathcal{E}(4,4) = 0\\
\end{array}
\]

%% G(x) = (1+erf(x/sqrt(2)))*(x/2) + exp(-x^2/2)/sqrt(2*pi)
%% etab = [[0]]
%% for n in range(1,21):
%%     row = [0]
%%     for k in range(1,n):
%%         row.append(N(etab[n-1][k-1] + G(etab[n-1][k]-etab[n-1][k-1]), digits=50))
%%     row.append(N(etab[n-1][n-1] + 0, digits=50))
%%     etab.append(row)

If $F(x) = x$ defines a uniform distribution on $[0,1]$
($\frac{dF}{dx} = 1$), then the expected value of $\max(X,t)$ is:
\[
G(t) = \frac{1}{2}(1+t^2)
\]
We get the following exact values:
\[
\begin{array}{c}
\mathcal{E}(0,0) = 0\\
\mathcal{E}(1,0) = 0\;\;\;,\;\;\; \mathcal{E}(1,1) = \frac{1}{2}\\
\mathcal{E}(2,0) = 0\;\;\;,\;\;\; \mathcal{E}(2,1) = \frac{5}{8}
\;\;\;,\;\;\; \mathcal{E}(2,2) = 1\\
\mathcal{E}(3,0) = 0\;\;\;,\;\;\; \mathcal{E}(3,1) = \frac{89}{128}
\;\;\;,\;\;\; \mathcal{E}(3,2) = \frac{153}{128}
\;\;\;,\;\;\; \mathcal{E}(3,3) = \frac{3}{2}\\
\mathcal{E}(4,0) = 0\;\;\;,\;\;\; \mathcal{E}(4,1) = \frac{24\,305}{32\,768}
\;\;\;,\;\;\; \mathcal{E}(4,2) = \frac{169}{128}
\;\;\;,\;\;\; \mathcal{E}(4,3) = \frac{57\,073}{32\,768}
\;\;\;,\;\;\; \mathcal{E}(4,4) = 2\\
\end{array}
\]

%% G(x) = (1+x^2)/2
%% etab = [[0]]
%% for n in range(1,21):
%%     row = [0]
%%     for k in range(1,n):
%%         row.append(etab[n-1][k-1] + G(etab[n-1][k]-etab[n-1][k-1]))
%%     row.append(etab[n-1][n-1] + 1/2)
%%     etab.append(row)

To give an example with a discrete distribution, if $X$ is defined by
rolling a fair $N$-sided die whose faces are numbered $1$ through $N$,
we have $F(x) = \frac{1}{N}\lfloor x\rfloor$ and $G(t) = \frac{N(N+1) -
  \lfloor t\rfloor(1+\lfloor t\rfloor-2t)}{2N}$
For $N=6$ faces, we get the following exact values:
\[
\begin{array}{c}
\mathcal{E}(0,0) = 0\\
\mathcal{E}(1,0) = 0\;\;\;,\;\;\; \mathcal{E}(1,1) = \frac{7}{2}\\
\mathcal{E}(2,0) = 0\;\;\;,\;\;\; \mathcal{E}(2,1) = \frac{17}{4}
\;\;\;,\;\;\; \mathcal{E}(2,2) = 7\\
\mathcal{E}(3,0) = 0\;\;\;,\;\;\; \mathcal{E}(3,1) = \frac{14}{3}
\;\;\;,\;\;\; \mathcal{E}(3,2) = \frac{49}{6}
\;\;\;,\;\;\; \mathcal{E}(3,3) = \frac{21}{2}\\
\mathcal{E}(4,0) = 0\;\;\;,\;\;\; \mathcal{E}(4,1) = \frac{89}{18}
\;\;\;,\;\;\; \mathcal{E}(4,2) = \frac{107}{12}
\;\;\;,\;\;\; \mathcal{E}(4,3) = \frac{215}{18}
\;\;\;,\;\;\; \mathcal{E}(4,4) = 14\\
\end{array}
\]

%% G(x) = (42-floor(x)*(1+floor(x)-2*x))/12
%% etab = [[0]]
%% for n in range(1,21):
%%     row = [0]
%%     for k in range(1,n):
%%         row.append(etab[n-1][k-1] + G(etab[n-1][k]-etab[n-1][k-1]))
%%     row.append(etab[n-1][n-1] + 7/2)
%%     etab.append(row)

For example, if we throw $n=10$ times an $N=6$-sided die and must
decide on the spot to keep-or-discard $k=5$ throws, we get a best
expected score of $\mathcal{E}(10,5) = \frac{181\,663}{7\,776} \approx
23.36201$ for their sum.

\hbox to\hsize{\hfill\hbox to8cm{\hrulefill}\hfill}

In contrast, let us now consider the following different situation,
still with $F$ being a fixed probability distribution on $\mathbb{R}$.
We draw $n$ i.i.d. values $X_0,\ldots,X_{n-1}$ from $F$ and we simply
take the $k$ best (i.e., we allow ourselves to wait until all values
are known to pick them).  Call $\mathcal{B}(n,k)$ the sum of these
$k$-best-of-$n$ values.  We wish to compute $\mathcal{B}(n,k)$.

Calling $Y_0\geq \cdots \geq Y_{n-1}$ the same values as
$X_0,\ldots,X_{n-1}$ but in decreasing order, the probability
$F_{n,\ell}(y)$ of the event $Y_\ell \leq y < Y_{\ell-1}$ is the
probability that exactly $\ell$ of the $X_i$ are $>y$ and exactly
$n-\ell$ are $\leq y$, namely:
\[
F_{n,\ell}(y) = \frac{n!}{\ell!\,(n-\ell)!} F(y)^{n-\ell} (1-F(y))^\ell
\]
so $\mathbb{P}(Y_j \leq y) = \sum_{\ell=0}^j F_{n,\ell}(y)$, and then
(\textit{“first computation”}):
\[
\begin{aligned}
\mathcal{B}(n,k) &= \sum_{j=0}^{k-1} \mathbb{E}(Y_j)
= \sum_{j=0}^{k-1} \sum_{\ell=0}^j \int y\,dF_{n,\ell}(y)
= \sum_{\ell=0}^{k-1} (k-\ell) \int y\,dF_{n,\ell}(y)\\
&= \sum_{\ell=0}^{k-1} \frac{n!(k-\ell)}{\ell!\,(n-\ell)!} \int y\,\Big((n-\ell)\,F(y)^{n-\ell-1} (1-F(y))^\ell - \ell F(y)^{n-\ell} (1-F(y))^{\ell-1}\Big)\,dF(y)\\
% &= \sum_{\ell=0}^{k-1} \frac{n!(k-\ell)}{\ell!\,(n-\ell)!} \int y\,\Big(\frac{n-\ell}{F(y)} - \frac{\ell}{1-F(y)}\Big)\,F(y)^{n-\ell} (1-F(y))^{\ell}\, dF(y)\\
&= \sum_{\ell=0}^{k-1} \frac{n!(k-\ell)}{\ell!\,(n-\ell-1)!} \int y\,F(y)^{n-\ell-1} (1-F(y))^{\ell}\, dF(y)\\
&\hphantom{=}\quad - \; \sum_{\ell'=0}^{k-2} \frac{n!(k-\ell'-1)}{\ell'!\,(n-\ell'-1)!} \int y\,F(y)^{n-\ell'-1} (1-F(y))^{\ell'}\, dF(y)\;\text{~[with $\ell'=\ell-1$]}\\
&= \sum_{\ell=0}^{k-1} \frac{n!}{\ell!\,(n-\ell-1)!} \int y\,F(y)^{n-\ell-1} (1-F(y))^{\ell}\, dF(y) = n \sum_{\ell=0}^{k-1} \mathcal{I}(n-1,\ell)\\
\end{aligned}
\]
where $\mathcal{I}(m,\ell) := \int y\,F_{m,\ell}(y)\, dF(y)$.

Alternatively, we may prefer to integrate by parts and write
(\textit{“second computation”}):
\[
\begin{aligned}
\mathcal{B}(n,k) &= \sum_{j=0}^{k-1} \mathbb{E}(Y_j)
= \sum_{j=0}^{k-1} \sum_{\ell=0}^j \int y\,dF_{n,\ell}(y)
= \sum_{\ell=0}^{k-1} (k-\ell) \int y\,dF_{n,\ell}(y)\\
&= \sum_{\ell=0}^{k-1} (k-\ell) \left([y F_{n,\ell}(y)]_{\mathrm{diff}}-\int F_{n,\ell}(y)\,dy\right)\\
\end{aligned}
\]
where the expression $[y F_{n,\ell}(y)]_{\mathrm{diff}}-\int
F_{n,\ell}(y)\,dy$ is to be read as $y_{\max} F_{n,\ell}(y_{\max}) -
y_{\min} F_{n,\ell}(y_{\min}) -\int_{y_{\min}}^{y_{\max}}
F_{n,\ell}(y)\,dy$ for some bounds $y_{\min},y_{\max}$ such that
$F_{n,\ell}$ is constant outside them (or as an appropriate limit of
this expression).

If $F(x) = \frac{1}{2}(1+\erf(\frac{x}{\sqrt{2}}))$ defines a standard
Gaussian distribution ($\frac{dF}{dx} =
\frac{1}{\sqrt{2\pi}}\exp(-\frac{x^2}{2})$), the computation of
$\mathcal{B}(n,k)$ in closed form is a delicate question in general,
nevertheless we can get (only expressing values for $2k\leq n$ since
we have $\mathcal{B}(n,n-k) = \mathcal{B}(n,k)$ here):
\[
\begin{array}{c}
\mathcal{B}(0,0) = 0\\
\mathcal{B}(1,0) = 0\;\;\;,\;\;\; \mathcal{B}(1,1) = 0\\
\mathcal{B}(2,0) = 0\;\;\;,\;\;\; \mathcal{B}(2,1) = \frac{1}{\sqrt{\pi}} \approx 0.56419\\
\mathcal{B}(3,0) = 0\;\;\;,\;\;\; \mathcal{B}(3,1) = \frac{3}{2\sqrt{\pi}} \approx 0.84628\\
\mathcal{B}(4,0) = 0\;\;\;,\;\;\; \mathcal{B}(4,1) = \frac{3}{\pi^{3/2}}(\pi-\arccos(\frac{1}{3})) \approx 1.02938
\;\;\;,\;\;\; \mathcal{B}(4,2) = \frac{6}{\pi^{3/2}}\arccos(\frac{1}{3}) \approx 1.32639\\
\end{array}
\]
(see \url{https://twitter.com/gro_tsen/status/1430529573171736579} and
appendix A in \url{https://rangevoting.org/BestVrange.html} for the
question of finding such closed form expressions, possibly involving
polylogarithms).

If $F(x) = x$ defines a uniform distribution on $[0,1]$
($\frac{dF}{dx} = 1$), then $\mathcal{I}(m,\ell) =
\frac{(m-\ell+1)}{(m+2)(m+1)}$ so $\mathcal{B}(n,k) =
\sum_{\ell=0}^{k-1} \frac{(n-\ell)(n-\ell-1)}{n+1} =
\frac{k(2n-k+1)}{2(n+1)}$.  In particular:
\[
\begin{array}{c}
\mathcal{B}(0,0) = 0\\
\mathcal{B}(1,0) = 0\;\;\;,\;\;\; \mathcal{B}(1,1) = \frac{1}{2}\\
\mathcal{B}(2,0) = 0\;\;\;,\;\;\; \mathcal{B}(2,1) = \frac{2}{3}
\;\;\;,\;\;\; \mathcal{B}(2,2) = 1\\
\mathcal{B}(3,0) = 0\;\;\;,\;\;\; \mathcal{B}(3,1) = \frac{3}{4}
\;\;\;,\;\;\; \mathcal{B}(3,2) = \frac{5}{4}
\;\;\;,\;\;\; \mathcal{B}(3,3) = \frac{3}{2}\\
\mathcal{B}(4,0) = 0\;\;\;,\;\;\; \mathcal{B}(4,1) = \frac{4}{5}
\;\;\;,\;\;\; \mathcal{B}(4,2) = \frac{7}{5}
\;\;\;,\;\;\; \mathcal{B}(4,3) = \frac{9}{5}
\;\;\;,\;\;\; \mathcal{B}(4,4) = 2\\
\end{array}
\]

For a discrete random variable, the “first computation” above fails to
be meaningful, but we can still perform the “second computation”
replacing the integral by the obvious discrete sum:
\[
\mathcal{B}(n,k)
= \sum_{\ell=0}^{k-1} (k-\ell) \left(y_{\max} F_{n,\ell}(y_{\max})
- y_{\min} F_{n,\ell}(y_{\min})
-\sum_{y=y_{\min}}^{y_{\max}} F_{n,\ell}(y)\right)
\]
(with $F(y_{\min})=0$ and $F(y_{\max})=1$ so that the two initial terms
are zero except when $\ell=0$ or $\ell=n$).

For an $N$-sided die with faces numbered from $1$ to $N$, this is
\[
\mathcal{B}(n,k) = \sum_{\ell=0}^{k-1} (k-\ell) \left(
(\text{$N$ if $\ell=0$ else $0$})
-\sum_{y=0}^{N-1} F_{n,\ell}(y)\right)
\]
This does not appear to have a nice closed form expression, but each
$\mathcal{B}(n,k)$ can be computed as a Laurent polynomial in $N$ :
\[
\begin{array}{c}
\mathcal{B}(0,0) = 0\\
\mathcal{B}(1,0) = 0\;\;\;,\;\;\; \mathcal{B}(1,1) = \frac{1}{2}(N+1)\\
\mathcal{B}(2,0) = 0\;\;\;,\;\;\; \mathcal{B}(2,1) = \frac{1}{6}(4N+3-N^{-1})
\;\;\;,\;\;\; \mathcal{B}(2,2) = N+1\\
\mathcal{B}(3,0) = 0\;\;\;,\;\;\; \mathcal{B}(3,1) = \frac{1}{4}(3N+2-N^{-1})\\
\;\;\;\ldots\;\;\; \mathcal{B}(3,2) = \frac{1}{4}(5N+4-N^{-1})
\;\;\;,\;\;\; \mathcal{B}(3,3) = \frac{3}{2}(N+1)\\
\mathcal{B}(4,0) = 0\;\;\;,\;\;\; \mathcal{B}(4,1) = \frac{1}{30}(24N+15-10N^{-1}+N^{-3})\\
\;\;\;\ldots\;\;\; \mathcal{B}(4,2) = \frac{1}{30}(42N+30-10N^{-1}-2N^{-3})\\
\;\;\;\ldots\;\;\; \mathcal{B}(4,3) = \frac{1}{30}(54N+45-10N^{-1}+N^{-3})
\;\;\;,\;\;\; \mathcal{B}(4,4) = 2(N+1)\\
\end{array}
\]
Of course the coefficient of the leading term in $N$ is the same value
tabulated above for the continuous uniform distribution.  For $N=6$
faces we get:
\[
\begin{array}{c}
\mathcal{B}(0,0) = 0\\
\mathcal{B}(1,0) = 0\;\;\;,\;\;\; \mathcal{B}(1,1) = \frac{7}{2}\\
\mathcal{B}(2,0) = 0\;\;\;,\;\;\; \mathcal{B}(2,1) = \frac{161}{36}
\;\;\;,\;\;\; \mathcal{B}(2,2) = 7\\
\mathcal{B}(3,0) = 0\;\;\;,\;\;\; \mathcal{B}(3,1) = \frac{119}{24}
\;\;\;,\;\;\; \mathcal{B}(3,2) = \frac{203}{24}
\;\;\;,\;\;\; \mathcal{B}(3,3) = \frac{21}{2}\\
\mathcal{B}(4,0) = 0\;\;\;,\;\;\; \mathcal{B}(4,1) = \frac{6797}{1296}
\;\;\;,\;\;\; \mathcal{B}(4,2) = \frac{6055}{648}
\;\;\;,\;\;\; \mathcal{B}(4,3) = \frac{15869}{1296}
\;\;\;,\;\;\; \mathcal{B}(4,4) = 14\\
\end{array}
\]

For example, if we throw $n=10$ times an $N=6$-sided die and keep the
best $k=5$ throws, we get an expected score of $\mathcal{B}(10,5) =
\frac{731\,015\,215}{30\,233\,088} \approx 24.17931$ for their sum.

\end{document}

See https://twitter.com/gro_tsen/status/1603363831715188736 for context.