9 Reinforcement Learning

“Information is the resolution of uncertainty.” (Claude Shannon)

Thus far we have discussed making decisions under uncertainty (Chapter 4 and Chapter 5) and two modes of data collection: field experiments and observational studies. In a field experiment we control the data-generation process before the study begins, whereas in an observational study we have no such control and must work with whatever data are produced.

What if we can choose which data to collect while the experiment is running? This leads to sequential (or adaptive) experimental design, in which each new observation is selected on the basis of the data gathered so far, creating a feedback loop between data generation and decision-making.

A loop in which a decision on the next sample collects data from the system under study, and observing the system feeds the next decision. — Figure 9.1: Sequential experimental design

This framework supports a wide range of applications and is implemented through several key algorithms. In this chapter, we will consider the most widely used among them: (i) Multi-Armed Bandits, (ii) Q-Learning, and (iii) Bayesian Optimization. One of the first practical demonstrations of reinforcement learning (then called trial-and-error learning) was Claude Shannon’s 1950s mechanical mouse Theseus, which learned to find its way through a maze (Shannon 1951).

9.1 Multi-Armed Bandits

When the number of alternatives is large and the testing budget limited, a different approach to A/B testing, the multi-armed bandit (MAB), can be more sample-efficient. MABs allow us to balance exploration (trying new things) and exploitation (sticking with what works) simultaneously, minimizing regret (the difference between our choice and the optimal one). MABs require that the outcome of each experiment is available immediately or with a small delay to decide on the next experiment (Steven L. Scott 2015). When time is of the essence, and there is an opportunity cost associated with delaying the decision, MABs are a better choice than traditional A/B testing.

The mathematical framework, the foundational Bayesian solution, Thompson sampling, and a comparison with classical A/B testing are covered in Chapter 6 (see Section 6.5.5). Here we focus on practical aspects: when to end experiments, how to handle contextual information, and design considerations.

Formally, there are $K$ alternatives (arms) and each arm $a$ has a reward distribution $f_a(y\mid \theta)$ and a value $v_a(\theta)$, its expected reward. The goal is to find the arm with the highest expected reward and accumulate the highest total reward in doing so. We can only observe the reward after we select an arm $a$, and we assume that we know the form of the distribution $f_a(y\mid \theta)$, where $a$ is the arm index, $y$ is the reward, and $\theta$ is the set of unknown parameters to be learned. Here are a few examples:

In online advertisements, we have $K$ possible ads to be shown and the probability of a user clicking is given by the vector $\theta = (\theta_1,\ldots,\theta_K)$ of success probabilities for $K$ independent binomial models, with $v_a(\theta) = \theta_a$. The goal is to find the ad with the highest click-through rate.
In website design, we have two design variables: the color of a button (red or blue) and its pixel size (27 or 40). We introduce two dummy variables, $x_c$ for color and $x_s$ for size, and the probability of a user clicking is given by \[ \mathrm{logit} P(\text{click} \mid \theta) = \theta_0 + \theta_c x_c + \theta_s x_s + \theta_{cs}x_cx_s, \] with each of the four combinations $(x_c,x_s)$ an arm and $v_a(\theta) = P(\text{click} \mid \theta, x_c, x_s)$.

Variations include controlling for background variables (covariates not under our control, such as time of day or user location) or replacing the binary outcome with a continuous variable (time spent on the website, money spent) using linear regression or another appropriate generalized linear model.

The key challenge is that acting greedily, always choosing the arm with the highest estimated reward $\hat a = \arg\max_a v_a(\hat \theta)$, may cause us to miss better alternatives. We are not sure about our estimates $\hat \theta$ and need to explore other options. Thompson sampling (Thompson 1933), the oldest and most elegant Bayesian solution, addresses this by sampling from the posterior distribution of each arm’s success probability and selecting the arm with the highest sample. This naturally balances exploration and exploitation: arms with high uncertainty have a chance to produce high samples and get selected, while arms with high posterior means are selected frequently. For Bernoulli rewards the algorithm achieves regret of order $\log T$ over $T$ steps, with a constant depending on the gaps between the arms, which matches the Lai and Robbins lower bound (Agrawal and Goyal 2012; Kaufmann, Korda, and Munos 2012).

When to End Experiments

Steps 2 and 3 of the Thompson sampling algorithm in Section 6.5.5 can be replaced by calculating the probability of an arm $a$ being the best $w_{at}$ and then choosing the arm by sampling from the discrete distribution $w_{1t},\ldots,w_{Kt}$. The probability of an arm $a$ being the best is given by \[ w_{at} = P(a \text{ is optimal } \mid y^t) = \int P(a \text{ is optimal } \mid \theta) P(\theta \mid y^t) d\theta, \] where $y^t = (y_1,\ldots,y_t)$ is the history of observations up to time $t$. We can calculate the probabilities $w_{at}$ using Monte Carlo. We can sample $\theta^{(1)},\ldots,\theta^{(G)}$ from the posterior distribution $p(\theta \mid y^t)$ and calculate the probability as \[ w_{at}\approx \dfrac{1}{G}\sum_{g=1}^G I(a = \arg\max_i v_i(\theta^{(g)})), \] where $I(\cdot)$ is the indicator function. This is simply the proportion of times the arm was the best in the $G$ samples.

Although using a single draw from posterior $p(\theta \mid y^t)$ (as in the original algorithm) is equivalent to sampling proportional to $w_{at}$, the explicitly calculated $w_{at}$ yields a useful statistic that can be used to decide on when to end the experiment.

We will use the regret statistic to decide when to stop. Regret is the difference in values between the truly optimal arm and the arm that is apparently optimal at time $t$. Although we cannot know the regret (it is unobservable), we can compute samples from its posterior distribution. We simulate the posterior distribution of the regret by sampling $\theta^{(1)},\ldots,\theta^{(G)}$ from the posterior distribution $p(\theta \mid y^t)$ and calculating the regret as \[ r^{(g)} = \max_a v_a(\theta^{(g)}) - v_{a^*_t}(\theta^{(g)}), \] where the first term is the value of the best arm within draw $g$, which varies from draw to draw, and $a^*_t$ is the champion, the single arm that is optimal most often across all draws, \[ a^*_t = \arg\max_a w_{at}. \] The two need to be kept apart: the regret is positive exactly on those draws in which some arm beats the champion.

Often, it is convenient to measure the regret on the percent scale, and then we use \[ r^{(g)} \leftarrow r^{(g)}/v_{a^*_t}(\theta^{(g)}) \]

We can demonstrate with a small example. The function below generates samples $\theta^{(g)}$

bandit <- function(x, n, alpha = 1, beta = 1, ndraws = 5000) {
  K <- length(x) # number of bandits
  prob <- matrix(nrow = ndraws, ncol = K)
  no <- n - x
  for (a in 1:K) { # posterior draws for each arm
    prob[, a] <- rbeta(ndraws, x[a] + alpha, no[a] + beta)
  }
  prob
}

Say we have three arms with 20, 30, and 40 sessions that have generated 12, 20, and 30 conversions. We assume a uniform prior for each arm $\theta_i \sim \text{Beta}(1,1)$ and generate 6 samples from the posterior of $\theta \mid y^t$.

x <- c(12, 20, 30)
n <- c(20, 30, 40)
set.seed(17) # Kharlamov
prob <- bandit(x, n, ndraws = 6)

Six posterior draws of the conversion rates of the three arms.
	$\theta_1$	$\theta_2$	$\theta_3$
1	0.60	0.63	0.58
2	0.62	0.62	0.74
3	0.69	0.53	0.67
4	0.49	0.59	0.73
5	0.61	0.51	0.69
6	0.47	0.64	0.69

Now, we calculate the posterior probabilities $w_{at} = P(a \text{ is optimal } \mid y^t)$ for each of the three arms

wat <- table(factor(max.col(prob), levels = 1:3)) / 6

Posterior probability that each arm is the best, from six draws.
1	2	3
0.17	0.17	0.67

Thus far, the third arm is the most likely to be optimal, with probability 67%. Now, we calculate the regret for each of the six draws from the posterior of $\theta \mid y^t$.

regret <- (apply(prob, 1, max) - prob[, 3]) / prob[, 3]

Regret of each of the six draws.
1	2	3	4	5	6
0.09	0	0.03	0	0	0

We compute the value row by row by subtracting the element in column 3 from the largest element of that row (arm 3 is the champion, because it has the highest chance of being optimal), and dividing by that same column 3 entry. All rows but 1 and 3 are zero. In the first row, the value is $(0.63-0.58)/0.58 = 0.09$, because column 2 is 0.05 larger than column 3. If we keep going down each row, we get a distribution of values that we could plot in a histogram. We can generate one for a larger number of draws (10000), redrawn from the same posteriors rather than continuing the six above.

prob <- bandit(x, n, ndraws = 10000)
regret <- (apply(prob, 1, max) - prob[, 3]) / prob[, 3]

The histogram of the value remaining in an experiment (regret). The vertical line is the 95th percentile, or the potential value remaining.

wat <- table(factor(max.col(prob), levels = 1:3)) / 10000

Posterior probability that each arm is the best, from 10000 draws.
1	2	3
0.08	0.2	0.72

Arm 3 has a 72% probability of being the best arm, so the value of switching away from arm 3 is zero in 72% of the cases. The 95th percentile of the value distribution is the potential value remaining in the experiment, a relative gap rather than a conversion rate, which in this case works out to be about 17%.

quantile(regret, 0.95)
##  95% 
## 0.17

You interpret this number as “We’re still unsure about the CvR for arm 3, but whatever it is, one of the other arms might beat it by as much as 17%.”

Google Analytics (Steven L. Scott 2013), for example, “ends the experiment when there’s at least a 95% probability that the value remaining in the experiment is less than 1% of the champion’s conversion rate. That’s a 1% improvement, not a one percentage point improvement. So if the best arm has a conversion rate of 4%, then we end the experiment if the value remaining in the experiment is less than .04 percentage points of CvR.”

Our three-arm experiment is far from that threshold. Its potential value remaining is 17%, well above 1%, so under this rule the experiment would keep running rather than declare arm 3 the winner.

Contextual Bandits

Traditional multi-armed bandit models, like the binomial model, assume independent observations with fixed reward probabilities. This works well when rewards are consistent across different groups and times. However, for situations with diverse user bases or fluctuating activity patterns, such as international audiences or browsing behavior, this assumption can be misleading.

For instance, companies with a global web presence may experience temporal effects as markets in Asia, Europe, and the Americas become active at different times of the day. Additionally, user behavior can change based on the day of the week, with people engaging in different browsing patterns and purchase behaviors. For example, individuals may research expensive purchases during work hours but make actual purchases on weekends.

Consider an experiment with two arms, A and B. Arm A performs slightly better during the weekdays when users browse but don’t make purchases, while Arm B excels during the weekends when users are more likely to make purchases. High traffic volume might lead a binomial model to declare Arm A the winner before observing any weekend behavior. This risk exists regardless of whether the experiment is conducted as a bandit or a traditional experiment. Bandit experiments, however, are particularly vulnerable to this bias due to their typically shorter durations.

To mitigate the risk of being misled by distinct sub-populations, two methods can be employed. If the specific sub-populations are known in advance or if there is a proxy for them, such as geographically induced temporal patterns, the binomial model can be adapted to a logistic regression model. \[ \mathrm{logit} P(\text{click}_a \mid \theta, x) = \beta_{0a} + \beta_a^\top x, \] where $x$ describes the circumstances or context of the observation. The success probability for selecting arm $a$ under the context $x$ is represented as $P(\text{click}_a \mid \theta, x)$. Each arm $a$ has its own specific coefficient denoted as $\beta_{0a}$ with one arm’s coefficient set to zero as a reference point. Additionally, there is another set of coefficients represented as $\beta_a$ that are associated with the contextual data and are learned as part of the model. With one $\beta$ shared by all arms, the difference between the logits of two arms would not depend on $x$, and the ranking of the arms could never change with the context. The value function can then be \[ v_a(\theta, x) = \mathrm{logit}^{-1}(\beta_{0a} + \beta_a^\top x). \]

If we lack knowledge about the crucial contexts, one option is to make the assumption that contexts are generated randomly from a context distribution. This approach is often exemplified by the use of a hierarchical model like the beta-binomial model. \[\begin{align*} \theta_{at} &\sim \text{Beta}(\alpha_a,\delta_a)\\ \text{click}_a \mid \theta_{at} &\sim \text{Bernoulli}(\theta_{at}), \end{align*}\] where $\theta = \{\alpha_a,\delta_a ~:~ a = 1,\ldots,K \}$, with value function $v_a(\theta) = \alpha_a/(\alpha_a + \delta_a)$.

Summary of MAB Experimentation

The design phase begins with defining your arms by identifying the different options you want to evaluate, such as different website layouts, pricing strategies, or marketing campaigns. Next, choose a bandit algorithm that balances exploration and exploitation in various ways. Popular choices include Epsilon-greedy, Thompson Sampling, and Upper Confidence Bound (UCB). Then set your parameters by configuring the algorithm parameters based on your priorities and expected uncertainty. For example, a higher exploration rate encourages trying new arms earlier.

During the analysis phase, track rewards by defining and measuring the reward metric for each arm, such as clicks, conversions, or profit. Monitor performance by regularly analyzing the cumulative reward and arm selection probabilities to see which arms are performing well and how the allocation strategy is adapting. Use statistical tools like confidence intervals or Bayesian methods to compare performance between arms and assess the significance of findings. Make adaptive adjustments by modifying the experiment based on ongoing analysis. You might adjust algorithm parameters, stop arms with demonstrably poor performance, or introduce new arms.

Set a clear stopping criterion by deciding when to end the experiment based on a predetermined budget, time limit, or desired level of confidence in the results.

9.2 Bellman Principle of Optimality

“An optimum policy has the property that whatever the initial state and initial decision are, the remaining decision sequence must be optimum for the state resulting from the first decision.” (Richard Bellman)

To solve sequential decision problems such as the secretary problem below, we rely on a fundamental concept in dynamic programming.

Example 9.1 (Secretary Problem) The Secretary Problem, also known as the marriage problem or sultan’s dowry problem, is a classic problem in which one selects the best option from a sequence of candidates (Gilbert and Mosteller 1966; Ferguson 1989). It is often framed as hiring a secretary, but it applies to any scenario, such as choosing a house or a spouse, in which options are evaluated sequentially and each decision is final.

In this problem, you receive $T$ offers and must either accept or reject the offer “on the spot.” You cannot return to a previous offer once you have moved on to the next one. Offers are in random order and can be ranked against those previously seen. The aim is to maximize the probability of choosing the offer with the greatest rank. There is an optimal $r$ ($1 \le r < T$) to be determined such that we examine and reject the first $r$ offers. Then of the remaining $T - r$ offers, we choose the first one that is best seen to date.

A decision strategy involves setting a threshold such that the first candidate above this threshold is hired, and all candidates below the threshold are rejected. The optimal strategy, known as the 37% rule, suggests that one should reject the first $r=T/e$ candidates and then select the first candidate who is better than all those seen so far.

The reasoning behind the 37% rule is based on the idea of balancing exploration and exploitation. By rejecting the first $T/e$ candidates, you gain a sense of the quality of the candidates but avoid committing too early. After that point, you select the first candidate who is better than the best among the initial $r$ candidates.

The 37% rule provides a probabilistic guarantee of selecting the best candidate with a probability close to 1/e (approximately 37%) as $T$ becomes large.

We use Bellman’s principle of optimality (Bellman 1957).

The solution to the secretary problem can be found via dynamic programming. Given an agent with utility function $u(x,d)$, with current state $x$, and decision $d$. The law of motion of $x_t$ is given by $x_{t+1} \sim p(\cdot \mid x_t,d_t)$. Bellman principle of optimality states that the optimal policy is given by the following recursion \[ V(x_t) = \max_{d_t} \left \{ u(x_t,d_t) + \gamma \E{V(x_{t+1})} \right \} \] where $\gamma$ is the discount factor. The optimal policy is given by \[ d_t^* = \arg \max_{d_t} \left \{ u(x_t,d_t) + \gamma \E{V(x_{t+1})} \right \}. \]

Now, back to the secretary problem. Let $y^t = (y_1,\ldots,y_t)$ denote the relative ranks of the first $t$ offers. State $x_t=1$ if the $t$th candidate is the best seen so far and $x_t=0$ otherwise. The decision $d_t=1$ if the $t$th candidate is hired and $d_t=0$ otherwise. Only hiring the overall best counts, so the utility is $u(x_t,d_t) = d_t\, I(t\text{th candidate is best overall})$, and it is worth noting that this is not the same as $x_t d_t$, which would pay for hiring any candidate who merely leads at the time. Before writing the Bellman equation we need the probability that a candidate who is best so far is best overall. Since “best of $T$” implies “best of the first $t$”, Bayes’ rule gives \[ P(\text{best of }T\mid x_t=1) = \dfrac{P(\text{best of }T)}{P(x_t=1)} = \dfrac{1/T}{1/t} = \dfrac{t}{T}. \] So the expected utility of hiring at a record at time $t$ is $\E{u \mid x_t=1, d_t=1} = t/T$, and this is the payoff that appears in the Bellman equation below. The $t$th offer is the best seen so far and places no restriction on the relative ranks of the first $t-1$ offers. Therefore, \[ p(x_t=1,y^{t-1}) = p(x_t=1)p(y^{t-1}) \] because the offers arrive in random order. Hence, we have \[ p(x_t=1 \mid y^{t-1}) = p(x_t=1) = \dfrac{1}{t}. \]

Let $p^*(x_{t-1}=0)$ be the probability of ending with the overall best candidate under the optimal strategy, given that the first $t-1$ offers have been passed over. It does not depend on whether offer $t-1$ was a record. Now we have to select the best candidate, given we have seen $t-1$ offers so far. The probability satisfies the Bellman equation \[ p^*(x_{t-1}=0) = \frac{t-1}{t} p^*(x_{t}=0) + \frac{1}{t}\max\left(t/T, p^*(x_{t}=0)\right). \] Once $t/T \ge p^*(x_t=0)$, that is, past the optimal threshold, the maximum is attained by $t/T$ and the recursion telescopes to \[ p^*(x_{t-1}=0) = \frac{t-1}{T} \sum_{\tau=t-1}^{T-1}\dfrac{1}{\tau}. \] Below the threshold it is optimal to keep looking, the maximum is attained by $p^*(x_t=0)$, and the value stays flat at whatever it reaches at the threshold. The expression above is therefore the value of the fixed policy that starts accepting at $t-1$, which is exactly how it is used next.

Remember, the strategy is to reject the first $r$ candidates and then select the first candidate better than all of them. The probability of selecting the best candidate is given by \[ P(\text{success}) = \dfrac{1}{T}\sum_{a=r+1}^T \dfrac{r}{a-1} = \dfrac{r}{T}\sum_{b=r}^{T-1} \dfrac{1}{b} \approx \dfrac{1}{T}\int_{r}^{T}\dfrac{r}{a} da = \dfrac{r}{T} \log \left ( \dfrac{T}{r} \right ). \] The summand is $r/(a-1)$, not $r/a$: the overall best sits at position $a$ with probability $1/T$, and it is then hired exactly when the best of the first $a-1$ fell in the rejected block, which has probability $r/(a-1)$. This is the same expression as the value function above with $t-1=r$. We optimize over $r$ by setting the derivative to zero: \[ \frac{\log \left(\frac{T}{r}\right)}{T}-\frac{1}{T} = 0, \] which gives the optimal $r=T/e$.

If we plug in $r=T/e$ back to the probability of success, we get \[ P(\text{success}) \approx \dfrac{1}{e} \log \left ( e \right ) = \dfrac{1}{e}. \]

Monte Carlo Simulations We start by showing how the secretary problem can be analyzed using simulations rather than the analytical derivations provided above.

# one trial: screen the first r - 1 candidates, then take the first one better than all of them
secretary <- function(n, r) {
  x <- sample(1:n, n)
  best_screen <- max(x[1:(r - 1)])
  later <- x[r:n]
  better <- which(later > best_screen)
  choice <- if (length(better) == 0) x[n] else later[min(better)]
  c(best = choice == max(x), quality = choice)
}

The sweep over the number of candidates screened is slow, so we ran it once and saved the table in data/secretary.csv, which the plots below read.

Figure 9.2: Probability of picking the best candidate (left) and average quality of the candidate chosen (right) against the number of candidates screened, out of 1000.

Figure 9.3: Probability of picking the best candidate (left) and average quality of the candidate chosen (right) against the number of candidates screened, out of 1000.

The left panel confirms the theoretical result: the probability of selecting the best candidate peaks close to the theoretical 37% near 370 screened out of 1000 candidates ($T/e \approx 368$), with Monte Carlo noise of about 1.5 points. The right panel reveals an interesting trade-off: average quality peaks much earlier, at roughly 20 to 60 candidates screened, then steadily declines. This happens because being too selective (screening too many) increases the risk of rejecting all strong candidates and being forced to accept the last one, regardless of quality. The two plots together illustrate the tension between maximizing the chance of finding the absolute best (left) and maximizing expected quality (right), objectives that lead to different optimal stopping rules.

9.3 Markov Decision Processes

A Markov Decision Process (MDP) models decision making in discrete time when outcomes are partly random and partly under the control of a decision maker. Almost all dynamic programming and reinforcement learning problems are formulated using the formalism of MDPs. Much of the early research on MDPs followed Ronald Howard’s 1960 book Dynamic Programming and Markov Processes (Howard 1960). In fact, the multi-armed bandit problem considered before is a special case of MDP with one state.

An MDP is defined by:

States ($S$): A set of states representing different scenarios or configurations. A key assumption is the Markov property: the future depends only on the current state and action, not on the history.
Actions ($A$): A set of actions available in each state.
Transition Probability ($P$): $P(s', r \mid s, a)$ is the probability of transitioning to state $s'$, receiving reward $r$, given that action $a$ is taken in state $s$.
Reward ($R$): A reward function $R(s, a, s')$ that gives the feedback signal immediately after transitioning from state $s$ to state $s'$, due to action $a$.
Discount Factor ($\gamma$): A factor between 0 and 1, which reduces the value of future rewards.

Mathematical Representation

The states $s_t$ and rewards $R_t$ in MDP are indexed by time $t$. The state at time $t+1$ is distributed according to the transition probability \[ P(s_{t+1}\mid s_t,a_t). \] The reward function is $R_s^a = \E{R_{t+1} \mid s, a}$.

The Markov property of the state is that the transition probability depends only on the current state and action and not on the history of states and actions. \[ P(s_{t+1}\mid s_t,a_t) = P(s_{t+1}\mid s_t,a_t,s_{t-1},a_{t-1},\ldots,s_0,a_0). \] The state is a sufficient statistic for the future.

In the case when the number of states is finite, we can represent the transition probability as a matrix $P_{ss'}^a = P(s_{t+1} = s' \mid s_t = s, a_t = a)$, where $s,s' \in S$ and $a \in A$. For a given action $a$, the transition probability matrix $P^a$ is a square matrix of size $|S| \times |S|$, where each row sums to 1 \[ P^a = \begin{bmatrix} P_{11}^a & P_{12}^a & \cdots & P_{1|S|}^a \\ P_{21}^a & P_{22}^a & \cdots & P_{2|S|}^a \\ \vdots & \vdots & \ddots & \vdots \\ P_{|S|1}^a & P_{|S|2}^a & \cdots & P_{|S||S|}^a \end{bmatrix} \] The reward function is also a matrix $R_s^a = \E{R_{t+1} \mid s_t = s, a_t = a}$.

Markov Reward Process

We can consider a simpler example, a Markov reward process. This is a special case of MDP when there is no action and the transition probability is simply a matrix $P_{ss'} = P(s_{t+1} = s' \mid s_t = s)$, where $s,s' \in S$. With no action to condition on, the transition probability matrix $P$ is a single square matrix of size $|S| \times |S|$ whose rows sum to 1.

Example 9.2 (Student Example) The graph below represents possible states (nodes) and transitions (links). Each node has a reward assigned to it, which corresponds to the reward function $R(s)$. The transition probabilities are shown on the links. The graph represents a Markov Chain where transitions are probabilistic and not controlled by an agent. The chain is adapted from the student example in David Silver’s UCL lectures on reinforcement learning (Silver 2015).

A Markov chain over the states Facebook, Class 1, Class 2, Class 3, Pub, Pass, Fail and Sleep, with a reward for each state and a transition probability on each arrow. — Figure 9.4: Student Example

There are no actions here, only transitions, so this is a Markov reward process rather than a full MDP. If we pick an initial state and sample a trajectory (a path on the graph above) by following the transition probabilities out of each state, we get a random walk on the graph. Sleep absorbs the walk: once there, it stays there and collects reward 0 forever. The reward for each state is shown in the graph. Counting from the state we currently occupy, the discounted value of the trajectory is \[ G_t = R(s_t) + \gamma R(s_{t+1}) + \gamma^2 R(s_{t+2}) + \cdots = \sum_{k=0}^\infty \gamma^k R(s_{t+k}), \] where $\gamma$ is the discount factor. The discount factor is a number between 0 and 1 that determines the present value of future rewards. A discount factor of 0 makes the agent myopic and only concerned about immediate rewards. A discount factor of 1 makes the agent strive for a long-term high reward. The discount of less than 1 is used to avoid infinite returns in cyclic Markov chains and allows us to discount less certain future rewards. The value of $\gamma$ is usually close to 1, for example 0.9 or 0.99. The value of $\gamma$ can be interpreted as the probability of the agent surviving from one time step to the next. We use the deliberately small $\gamma = 0.5$ below so that the arithmetic stays legible and short trajectories capture most of the return.

We can calculate sample returns $G_t$ for this Markov Chain. We first read in the reward vector

Rewards
	Facebook	Class 1	Class 2	Class 3	Pub	Pass	Fail	Sleep
Reward	-1	-2	-2	-2	3	10	-20	0

and then the transition probability matrix

Transition probability matrix of the Student example.
	Facebook	Class 1	Class 2	Class 3	Pub	Pass	Fail	Sleep
Facebook	0.9	0.1
Class 1	0.5		0.5
Class 2				0.8				0.2
Class 3					0.4	0.6
Pub		0.2	0.3	0.3			0.2
Pass								1
Fail								1
Sleep								1

Now we check that all of the rows sum to 1

Transition probability matrix row sums.
Facebook	Class 1	Class 2	Class 3	Pub	Pass	Fail	Sleep
1	1	1	1	1	1	1	1

Given the transition probability matrix p, we can sample possible trajectories. First, we define a tr(s,m) convenience function that generates a trajectory of length m starting from state s

set.seed(17)
# Sample column s' using probabilities in row s
jump <- function(s) sample(rownames(p), 1, prob = p[s, ])
# Function to generate a trajectory
tr <- function(s, m) {
  traj <- c(s)
  for (i in 1:m) {
    traj <- c(traj, jump(traj[i]))
  }
  return(traj)
}

Now, we generate 6 trajectories of length 5 starting from state “Pub”

Six sample trajectories of length 5 from Pub.
Pub	Class 3	Pub	Class 2	Class 3	Pass
Pub	Class 2	Class 3	Pass	Sleep	Sleep
Pub	Class 2	Class 3	Pub	Fail	Sleep
Pub	Fail	Sleep	Sleep	Sleep	Sleep
Pub	Fail	Sleep	Sleep	Sleep	Sleep
Pub	Class 3	Pass	Sleep	Sleep	Sleep

Now we can calculate the discounted value $G_t$ of each of the trajectories

Discounted value of each trajectory.
2.7	2.8	0.62	-7	-7	4.5

We can calculate the discounted value for 1000 trajectories, where getR(s) returns the reward of state s from the reward table above

# Value function of a trajectory
value <- function(s, m, gamma = 0.5) {
  traj <- tr(s, m)
  disc <- gamma^(0:m)
  return(sum(sapply(traj, getR) * disc))
}
vpub <- replicate(1000, value("Pub", 6))
hist(vpub, main = "", xlab = "Discounted reward")

Discounted value of 1000 trajectories from Pub.

mean(vpub)
## 0.72

We can see that the discounted rewards fall into two clusters, and what separates them is not whether a path reaches “Fail” but when. The left cluster is the single value $-7$, which has probability $0.2$: these are the paths that go from Pub straight to Fail on the first step, each worth $3 + 0.5 \times (-20) = -7$. All other paths are worth between $0.25$ and $4.5$, including the 5.7% that reach Fail later, because by then the $-20$ is discounted by a factor of $0.125$ or less. This right cluster has two bumps: a spike at $4.5$ with probability $0.18$, the paths that go from Pub to Class 3 to Pass and then to Sleep, worth $3 + 0.5 \times (-2) + 0.25 \times 10 = 4.5$, and a broad bump between $1$ and $3$ that holds 55% of the paths.

This average is a Monte Carlo estimate of a return truncated at seven states, not the infinite-horizon value. We can get the value exactly. Writing the Bellman equation for every state at once as $V = R + \gamma P V$ and solving, \[ V = (I - \gamma P)^{-1} R, \]

gamma <- 0.5
V <- solve(diag(nrow(p)) - gamma * p, R[rownames(p), ])
names(V) <- rownames(p)
knitr::kable(t(V), digits = 2)

Exact state values.
Facebook	Class 1	Class 2	Class 3	Pub	Pass	Fail	Sleep
-2.1	-2.9	-1.6	1.1	0.65	10	-20	0

Truncating at seven states costs almost nothing here: the exact seven-state return from Pub is $0.648$ against an exact value of $0.646$, a bias of about $0.001$, because at $\gamma = 0.5$ the tail is worth very little even for the 22% of paths still short of the absorbing Sleep state. The remaining distance to the simulated $0.72$ is Monte Carlo error, not truncation: the standard error of an average over 1000 paths is about $0.12$ here.

The value of a state is the expected discounted reward starting from that state \[ V(s) = \E{G_t \mid s_t = s}. \] It evaluates the long-term value of state $s$ (the goodness of a state). It can be drastically different from the reward value associated with the state. In our student example, the reward for the “Pub” state is 3, but the value is only 0.65, because a fifth of the time the Pub leads to Fail and its reward of $-20$.

The value of a state can be calculated recursively using the Bellman equation \[\begin{align*} V(s) &= \E{G_t \mid s_t = s} \\ &= \E{R(s_t) + \gamma R(s_{t+1}) + \gamma^2 R(s_{t+2}) + \cdots \mid s_t = s} \\ &= \E{R(s_t) + \gamma G_{t+1} \mid s_t = s} \\ &= R(s) + \gamma\,\E{V(s_{t+1}) \mid s_t = s} \\ &= \sum_{s'} P(s' \mid s) \left[R(s) + \gamma V(s')\right], \end{align*}\] where the last step uses $\sum_{s'} P(s' \mid s) = 1$ to pull the known $R(s)$ back inside the sum. In matrix form this is $V = R + \gamma P V$, the system solved in Example 9.2.

Example 9.3 (MDP for Forest Management) We can consider one of the classic examples of a Markov Decision Process (MDP). Imagine you need to calculate an optimal policy to manage a forest to prevent possible fires. The goal is to decide between two possible actions to either ‘Wait’ or ‘Cut’. They correspond to balancing between ecological preservation and economic gain, considering the random event of a fire. We can break down the elements of this model.

States: Represent the age of the forest. The states are denoted as $\{1, 2,\ldots, S\}$, where 1 is the youngest state (just after a fire or cutting), and $S$ is the oldest state of the forest.
Actions: There are two actions available:
- ‘Wait’ (Action 1): Do nothing and let the forest grow for another year.
- ‘Cut’ (Action 2): Harvest the forest, which brings immediate economic benefit but resets its state to the youngest.
Probabilities: There’s a probability ‘p’ each year that a fire occurs, regardless of the action taken. If a fire occurs, the forest returns to state 1.
Transition Matrix (P): This matrix defines the probabilities of moving from one state to another, given a specific action.

We build the transition probability array and the reward matrix of the Forest example in R, for four states, a fire probability of $p = 0.01$ per year, a reward $r_1 = 10$ for waiting in the oldest state and a reward $r_2 = 1$ for harvesting once the forest has grown beyond the youngest state.

fire <- 0.01
P_forest <- array(0, c(4, 4, 2))
P_forest[cbind(1:4, c(2, 3, 4, 4), 1)] <- 1 - fire  # Wait: the forest grows a year older
P_forest[, 1, 1] <- P_forest[, 1, 1] + fire         # a fire sends it back to state 1
P_forest[, 1, 2] <- 1
R_forest <- cbind(Wait = c(0, 0, 0, 10), Cut = c(0, 1, 1, 1))

The array $P$ has size $(|S| \times |S| \times |A|)$, so that P_forest[s, s', a] is $P(s' \mid s, a)$ and the two action slices are P_forest[,,1] and P_forest[,,2]. There are four states, $S = \{1,2,3,4\}$, and two actions $A = \{1,2\}$. The transition probability matrices for each action are: \[ P(\cdot \mid s, a=1) = \begin{bmatrix} 0.01 & 0.99 & 0.00 & 0.00 \\ 0.01 & 0.00 & 0.99 & 0.00 \\ 0.01 & 0.00 & 0.00 & 0.99 \\ 0.01 & 0.00 & 0.00 & 0.99 \end{bmatrix}, \quad P(\cdot \mid s, a=2) = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \end{bmatrix} \]

The reward matrix $R$ of size $|S| \times |A|$ specifies the immediate reward for taking each action in each state:

\[ R = \begin{bmatrix} 0 & 0 \\ 0 & 1 \\ 0 & 1 \\ 10 & 1 \end{bmatrix} \]

With these transition probabilities and rewards defined, we can solve for the optimal policy using either value iteration or policy iteration. Value iteration works by iteratively updating the value function until convergence, while policy iteration alternates between policy evaluation and policy improvement steps. In exact arithmetic both algorithms converge to the same optimal policy, though they may differ in computational efficiency depending on the problem structure.

Value iteration takes a few lines of R. Starting from $V = 0$, it applies the update $V(s) \leftarrow \max_a \left\{ R(s,a) + \gamma \sum_{s'} P(s' \mid s, a) V(s') \right\}$ until the largest change in $V$ is below a tolerance, and then returns the value function and the action that attains the maximum in each state.

value_iteration <- function(P, R, gamma, tol = 1e-10) {
  V <- rep(0, nrow(R))
  repeat {
    Q <- sapply(1:ncol(R), function(a) R[, a] + gamma * P[, , a] %*% V)
    V_new <- apply(Q, 1, max)
    if (max(abs(V_new - V)) < tol) break
    V <- V_new
  }
  list(V = V_new, policy = max.col(Q, ties.method = "first"))
}
vi <- value_iteration(P_forest, R_forest, gamma = 0.9)
cat("Optimal policy:", vi$policy, "\nValues:", sprintf("%.1f", vi$V), "\n")
## Optimal policy: 1 1 1 1 
## Values: 70.7 78.7 87.6 97.6

For this simple MDP with a discount factor of $\gamma = 0.9$, the optimal policy is $\pi^* = (1, 1, 1, 1)$: the Wait action is optimal in every state, and the values of the four states are 70.7, 78.7, 87.6 and 97.6. Because a fire strikes with probability only $p = 0.01$ while the mature-forest reward $r_1 = 10$ dwarfs the harvest reward $r_2 = 1$, cutting is never worthwhile; letting the forest keep growing and collecting $r_1$ repeatedly in the oldest state yields a higher expected discounted reward than harvesting.

Example 9.4 (Game of Chess as an MDP) We can view chess as an MDP.

In chess, a state $s$ represents the configuration of the chessboard at any given time. This includes the positions of all the pieces (pawns, knights, bishops, rooks, queen, and king) for both players (white and black). The arrangement alone is not a Markov state: the identical placement admits castling in one game and not in another, depending on whether the king or rook has already moved. A state must therefore record the placement together with the side to move, the remaining castling rights, the en-passant square if any, the halfmove counter for the fifty-move rule, and the positions that have occurred since the last irreversible move, with their counts (a draw may be declared when the same position recurs). With those components the state does determine the legal moves and the game is Markov. The game starts in a standard initial state (the standard chess setup) and progresses through a series of states as moves are made. If the game is played to completion, it ends in a terminal state (checkmate or draw). In a timed game, the game can also end when a player runs out of time.

Actions $a$ in chess are the legal moves that can be made by the player whose turn it is to move. This includes moving pieces according to their allowed movements, capturing opponent pieces, and special moves like castling or en passant. The set of actions available changes with each state, depending on the position of the pieces on the board.

In chess, the transition probability is deterministic for the most part, meaning that the outcome of a specific action (move) is certain and leads to a predictable next state. For example, moving a knight from one position to another (assuming it’s a legal move) will always result in the same new configuration of the chessboard. However, in the context of playing against an opponent, there is uncertainty in predicting the opponent’s response, which can be seen as introducing a probabilistic element in the larger view of the game. The game is an MDP only against an opponent who follows a fixed policy, which then enters the transition probabilities.

Defining a reward function $R$ in chess can be complex. In the simplest form, the reward could be associated with the game’s outcome: a win, loss, or draw. Wins could have positive rewards, losses negative, and draws could be neutral or have a small positive/negative value. Alternatively, more sophisticated reward functions can be designed to encourage certain strategies or positions, like controlling the center of the board, protecting the king, or capturing opponent pieces.

Chess is a game of perfect information, meaning all information about the game state is always available to both players. While the number of states in chess is finite, it is extremely large, making exhaustive state analysis (like traditional MDP methods) computationally impractical.

In practice, solving chess as an MDP, especially using traditional methods like value iteration or policy iteration, is not feasible due to the enormous state space. Modern approaches involve heuristic methods, machine learning, and deep learning techniques. For instance, advanced chess engines and AI systems like AlphaZero (Silver et al. 2017) use deep neural networks and reinforcement learning to evaluate board positions and determine optimal moves, but they do not solve the MDP in the classical sense.

The goal in an MDP is to find a policy $\pi(a \mid s)$, a distribution over actions in each state, of which a deterministic policy $a = \pi(s)$ is the special case putting all mass on one action, that maximizes the sum of discounted rewards: \[ V^\pi(s) = \E[\pi]{G_t \mid S_t = s}, \] where \[ G_t = \sum_{k=0}^{\infty} \gamma^k R(s_{t+k}, a_{t+k}, s_{t+k+1}) \] with actions drawn as $a_{t+k} \sim \pi(\cdot \mid s_{t+k})$.

Function $V^\pi(s)$ is the value of state s under policy $\pi$. Similarly we can define the action-value function $Q^\pi(s,a)$ as the value of taking action $a$ in state $s$ under policy $\pi$: \[ Q^\pi(s,a) = \E[\pi]{G_t \mid S_t = s, A_t = a}. \]

Bellman Equations for MDP simply state that the value of a state is the sum of the immediate reward and the discounted value of the next state \[ V^\pi(s) = \E[\pi]{R_{t+1} + \gamma V^{\pi}(S_{t+1})\mid S_t = s} = \sum_{a\in A}\pi(a\mid s)\left(R_s^a + \gamma \sum_{s'\in S}P^a_{ss'}V^\pi(s') \right). \] The action-value function satisfies the following Bellman equation \[ Q^\pi(s,a) = \E[\pi]{R_{t+1} + \gamma Q^{\pi}(S_{t+1}, A_{t+1})\mid S_t = s, A_t = a}. \] The value function can be defined as an expectation over the action-value function \[ V^\pi(s) = \E[\pi]{Q^\pi(s,a)\mid S_t = s} = \sum_{a\in A}\pi(a\mid s)Q^\pi(s,a). \] Expanding one step, \[ Q^\pi(s,a) = R_s^a + \gamma \sum_{s'\in S}P_{ss'}^a V^\pi(s') = R_s^a + \gamma \sum_{s'\in S}P_{ss'}^a\sum_{a'\in A}\pi(a'\mid s')Q^\pi(s',a'). \] Now we can define the Bellman equation in the matrix form \[ V^\pi = R^\pi + \gamma P^\pi V^\pi, \] where $R^\pi(s) = \sum_a \pi(a\mid s) R_s^a$ and $P^\pi_{ss'} = \sum_a \pi(a\mid s) P^a_{ss'}$. The direct solution is then \[ V^\pi = (I - \gamma P^\pi)^{-1}R^\pi. \] The optimal value function $V^*(s)$ is the value function for the optimal policy $\pi^*(s)$ \[ V^*(s) = \max_\pi V^\pi(s). \] The optimal action-value function $Q^*(s,a)$ is the action-value function for the optimal policy $\pi^*(s)$ \[ Q^*(s,a) = \max_\pi Q^\pi(s,a). \] The optimal policy $\pi^*(s)$ is the policy that maximizes the value function \[ \pi^*(s) = \arg\max_a Q^*(s,a). \] The optimal value function satisfies the Bellman optimality equation \[ V^*(s) = \max_a Q^*(s,a). \] and vice versa \[ Q^*(s,a) = R_s^a + \gamma \sum_{s'\in S}P_{ss'}^a V^*(s'). \]

The Bellman optimality equation is non-linear and is typically solved iteratively using Value Iteration, Policy Iteration, or Q-learning, which we return to below.

Example 9.5 (Q-Values and Deal or No Deal) Deal or No Deal is a popular TV show (Post et al. 2008) where a contestant is presented with a number of sealed boxes, each containing a prize. The contestant selects a box and then proceeds to open the remaining boxes one by one. After a certain number of boxes have been opened, the banker makes an offer to buy the contestant’s box. The contestant can either accept the offer and sell the box or reject the offer and continue opening boxes. The game continues until the contestant either accepts an offer or opens all the boxes. The goal is to maximize the expected utility of the amount received, either the banker’s offer if the contestant deals, or the prize in their own box if they refuse every offer. The rule of thumb is to continue as long as there are two large prizes left, because the continuation value is large. For example, with three prizes and two large ones, risk averse people will naively choose deal, when if they incorporated the continuation value they would choose no deal.

Let $s$ denote the current state of the system and $a$ an action. The $Q$-value, $Q_t(s,a)$, is the value of using action $a$ today and then proceeding optimally in the future. We use $a=1$ to mean no deal, and $a=0$ means deal. The Bellman equation for $Q$-values becomes \[ Q_{t} ( s , a) = u( s , a ) + \sum_{ s^\star } P( s^\star \mid s ,a ) \max_{ a' } Q_{t+1} ( s^\star , a' ) \] where $u(s,a)$ is the immediate utility of taking action $a$ in state $s$. Because the transition model is fully specified here, these $Q$-values are computed exactly by backward induction, in contrast to the model-free Q-learning introduced later in the chapter, which estimates them from sampled episodes. The value function and optimal action are given by \[ V(s) = \max_a Q ( s , a ) \; \; \text{and} \; \; a^\star = \arg\max_a Q ( s , a ) \]

Transition Matrix: Consider the problem where you have three prizes left. Now $s$ is the current state of three prizes. \[ s^\star = \{ \text{all sets of two prizes} \} \; \; \text{and} \; \; P( s^\star \mid s, a =1) = \frac{1}{3} \] where the transition matrix is uniform to the next state. There’s no continuation for $P( s^\star \mid s, a =0)$.

Utility: The utility of the next state depends on the contestant’s value for money and the bidding function of the banker \[ u( B ( s^\star ) ) = \frac{ B ( s^\star )^{1-\eta} -1 }{1 - \eta } \] in the power utility case, where $\eta$ is the coefficient of relative risk aversion. This is the same parameter written $\gamma$ in Example 4.3; we rename it here only to avoid colliding with the discount factor $\gamma$ used throughout this chapter.

Expected value implies $B( s ) = \bar{s}$ where $s$ are the remaining prizes.

We assume the following bidding function for the banker. With three prizes left: \[ B( s) = 0.305 \times \text{big} + 0.5 \times \text{small} \] and with two prizes left \[ B( s) = 0.355 \times \text{big} + 0.5 \times \text{small} \]

Three prizes left: $s = \{ 750 , 500 , 25 \}$.

Assume the contestant is risk averse with log-utility $U(x) = \log x$, which is the $\eta \to 1$ limit of the power utility above. Risk aversion is what drives the example: under expected value alone the two branches below are exactly indifferent, both worth 425. If the banker offers the expected value, we get \[ u( B( s = \{ 750 , 500 , 25 \}) ) = \log ( 1275/3 ) = 6.052 \] and so $Q_t ( s , a= 0 ) = 6.052$.

In the continuation problem, $s^\star = \{ s_1^\star , s_2^\star , s_3^\star \}$ where $s_1^\star = \{750,500 \}$ and $s_2^\star = \{ 750,25 \}$ and $s_3^\star = \{ 500,25 \}$.

We’ll have offers $625 , 387.5 , 262.5$ under the expected value. As the banker offers expected value, the optimal action at time $t+1$ is to take the deal $a=0$ with Q-values given by \[\begin{align*} Q_{t} ( s , a=1) & = \sum_{ s^\star } P( s^\star \mid s ,a =1) \max_{ a } Q_{t+1} ( s^\star , a ) \\ & = \frac{1}{3} \left ( \log (625) + \log (387.5) + \log (262.5) \right ) = 5.989 \end{align*}\] as immediate utility $u( s,a ) = 0$. Hence as \[ Q_{t} ( s , a=1)=5.989 < 6.052 = Q_{t} ( s , a=0) \] the optimal action is $a^\star = 0$, deal. Continuation value is not large enough to overcome the generous (expected value) offered by the banker.

Sensitivity analysis: we perform it by assuming different Banker’s bidding function. If we use the two-prize bidding function: \[ B( s) = 0.355 \times \text{big} + 0.5 \times \text{small}, \] Hence \[\begin{align*} B( s_1^\star = \{750,500 \}) & = 516.25 \\ B( s_2^\star = \{ 750,25 \}) & = 278.75 \\ B( s_3^\star = \{ 500,25 \}) & = 190 \end{align*}\]

The optimal action with two prizes left for the contestant is \[\begin{align*} Q_{t+1} ( s_1^\star , a=1) & = \frac{1}{2} \left ( \log (750) + \log (500) \right ) = 6.417 \\ & > 6.246 = Q_{t+1} ( s_1^\star , a=0) = \log \left ( 516.25 \right ) \\ Q_{t+1} ( s_2^\star , a=1) & = \frac{1}{2} \left ( \log (750) + \log (25) \right ) = 4.9194 \\ & < 5.63 = Q_{t+1} ( s_2^\star , a=0) = \log \left ( 278.75 \right ) \\ Q_{t+1} ( s_3^\star , a=1) & = \frac{1}{2} \left ( \log (500) + \log (25) \right ) = 4.717 \\ & < 5.247 = Q_{t+1} ( s_3^\star , a=0) = \log \left ( 190 \right ) \\ \end{align*}\] Hence the future optimal policy will be no deal under $s_1^\star$, and deal under $s_2^\star , s_3^\star$.

Therefore solving for $Q$-values at the previous step gives \[\begin{align*} Q_{t} ( s , a=1) & = \sum_{ s^\star } P( s^\star \mid s ,a =1) \max_{ a } Q_{t+1} ( s^\star , a ) \\ & = \frac{1}{3} \left ( 6.417+ 5.63 + 5.247 \right ) = 5.765 \end{align*}\] with a monetary equivalent as $\exp(5.765 ) = 318.9$.

With three prizes, we have \[\begin{align*} Q_{t} ( s , a=0) & = u( B( s = \{ 750 , 500 , 25 \}) ) \\ & = \log \left ( 0.305 \times 750 + 0.5 \times 25 \right ) \\ & = \log ( 241.25 ) = 5.48. \end{align*}\] The contestant is offered $ 241.

Now we have $Q_{t} ( s , a=1)= 5.765 > 5.48 = Q_{t} ( s , a=0)$ and the optimal action is $a^\star = 1$, no deal. The continuation value is large. The offer is $241 against a certainty equivalent of $319, a 32% premium for continuing.

MDP Solvers

The underlying approach behind all MDP solvers is to iteratively apply the Bellman equations until convergence. The main difference between the solvers is how they update the value function. All of them use a dynamic programming approach to find an optimal policy. Dynamic programming is a method for solving complex problems by breaking them down into simpler subproblems. It is applicable to problems exhibiting the properties of overlapping subproblems and optimal substructure. If a problem can be solved by combining optimal solutions to non-overlapping subproblems, the strategy is called divide and conquer instead. The Bellman equation gives the recursion, and the value function stores the solutions of the subproblems, which is why dynamic programming applies to MDPs.

First, we consider how to find the values of states under a given policy $\pi$. We can iteratively apply Bellman expectation backup. We update the values using the following update rule \[ V_{k+1}(s) = \sum_{a} \pi(a \mid s) \sum_{s'} P(s' \mid s, a)[R(s, a, s') + \gamma V_k(s')]. \] We will introduce the short-cut notation \[ P_{ss'}^a = P(s' \mid s, a), \quad R_s^a = \sum_{s'} P(s' \mid s, a)R(s, a, s'). \] Then in matrix form the update rule becomes \[ V_{k+1} = R^{\pi} + \gamma P^{\pi} V_k. \]

Policy Iteration

The policy iteration algorithm involves two main steps: policy evaluation and policy improvement, which are iteratively applied until convergence. We start with an arbitrary policy $\pi_0$. \[\begin{align*} V^{\pi_k} &= R^{\pi_k} + \gamma P^{\pi_k} V^{\pi_k} = (I - \gamma P^{\pi_k})^{-1} R^{\pi_k} \\ \pi_{k+1}(s) &= \arg\max_a \left( R^a + \gamma P^a V^{\pi_k} \right) = \arg\max_a Q^{\pi_k}(s,a) \end{align*}\] The evaluation step is run to convergence (or solved directly, as the inverse above shows) before each improvement; stopping the evaluation after a single sweep gives a different algorithm, modified policy iteration, with a different iteration count. The last step is to simply choose the action that maximizes the expected return in each state. Convergence is guaranteed: each improvement step yields a policy at least as good (the policy improvement theorem), and a finite MDP has finitely many deterministic policies, so policy iteration terminates after finitely many improvements. We stop when a full improvement pass leaves the policy unchanged.

The Bellman optimality equation is a fundamental part of finding the best policy in MDPs. \[ V^*(s) = \max_a \sum_{s', r} P(s', r \mid s, a)[r + \gamma V^*(s')] \] The optimal policy is then \[ \pi^*(s) = \arg\max_a \sum_{s', r} P(s', r \mid s, a)[r + \gamma V^*(s')] \]

Given an optimal policy, we can subdivide it into two parts: the optimal first action $A^*$ and the optimal policy from the next state $S'$. The optimal value $V^*$ can be found using one-step lookahead \[ V^*(s) = \max_a \left\{ R_s^a + \gamma \sum_{s'\in S} P_{ss'}^a V^*(s') \right\} \]

Value Iteration

This allows us to define another approach to solving MDPs, called value iteration. The value iteration algorithm starts with an arbitrary value function and iteratively applies the Bellman optimality backup. The algorithm updates the value function using the following update rule \[ V_{k+1}(s) = \max_a \left\{ R_s^a + \gamma \sum_{s'\in S} P_{ss'}^a V_k(s') \right\}. \] In matrix form, the update rule becomes \[ V_{k+1} = \max_a \left\{ R^a + \gamma P^a V_k \right\}. \] The algorithm stops when the maximum change in the value function is below a threshold. The optimal policy can be found by maximizing the right-hand side of the Bellman optimality equation \[ \pi^*(s) = \arg\max_a \left\{ R_s^a + \gamma \sum_{s'\in S} P_{ss'}^a V^*(s') \right\}. \]

In practice, exactly solving the Bellman Expectation Equation in the policy evaluation step can be computationally expensive for large state spaces. Approximate methods may be used.

Example 9.6 (MDP for a Maze) We use a mazemdp archive by Sally Gao, Duncan Rule, Yi Hao to demonstrate the value and policy iterations. Both are applied to the problem of finding an optimal policy for a maze. The maze is represented as a grid, with each cell either being a wall or empty. Both solvers are model-based, so the agent (decision maker) knows the layout and the transition structure, and each computes an optimal action for every cell rather than a single path; following the greedy actions from any starting cell then traces an optimal route to the goal in the top left corner. The agent can move up, down, left, or right, but not diagonally (actions). Moving into a wall keeps the agent in the same cell. Every cell carries a reward of $-1$ except the goal, which carries $+10$, and rewards are discounted with $\gamma=0.9$. Discounting is what makes short paths preferable: a cell $d$ steps from the goal has value $-10 + 110(0.9)^d$, strictly decreasing in $d$. The goal is not absorbing in this implementation, since bumping into the outer wall keeps the agent in place, so $V(\text{goal}) = 10/(1-\gamma) = 100$. The goal is to find the optimal policy that maximizes the expected return (sum of discounted rewards) for the agent.

Figure 9.5 and Figure 9.6 show snapshots from policy and value iterations.

Policy iteration converged after 23 iterations, the sweep at which a full pass leaves the policy unchanged.

Model-Free Methods

Both policy and value iterations we’ve considered thus far assume that transition probabilities between states given actions are known. However, this is often not the case in many real-world problems. Model-free methods learn through trial and error, by interacting with the environment and observing the rewards (Sutton and Barto 2018). The first method we consider is Monte Carlo methods. Monte Carlo methods for Markov Decision Processes (MDPs) are a class of algorithms used for finding optimal policies when the model of the environment (i.e., the transition probabilities and rewards) is unknown or too complex to model explicitly. These methods rely on learning from experience, specifically from complete episodes of interaction with the environment.

Generate Episodes: An episode is a sequence of states, actions, and rewards, from the start state to a terminal state. \[ S_0, A_0, R_1, S_1, A_1, R_2, \ldots, S_{T-1}, A_{T-1}, R_T \sim \pi. \] In Monte Carlo methods, these episodes are generated through actual or simulated interaction with the environment, based on a certain policy.
Estimate Value Functions: Unlike dynamic programming methods, which update value estimates based on other estimated values, Monte Carlo methods update estimates based on actual returns received over complete episodes. This involves averaging the returns received after visits to each state. We use empirical mean to estimate the expected value.
Policy Improvement: After a sufficient number of episodes have been generated and value functions estimated, the policy is improved based on these value function estimates.

Monte Carlo methods require complete episodes to update value estimates. This makes them particularly suitable for episodic tasks, where interactions naturally break down into separate episodes with clear starting and ending points. MC methods require sufficient exploration of the state space. This can be achieved through various strategies, like $\epsilon$-greedy policies, where there’s a small chance of taking a random action instead of the current best-known action. In this case, the policy is given by \[ \pi(a \mid s) = \begin{cases} 1 - \epsilon + \frac{\epsilon}{|A|} & \text{if } a = \arg\max_{a'} Q(s,a') \\ \frac{\epsilon}{|A|} & \text{otherwise} \end{cases} \] where $\epsilon$ is the probability of taking a random action and $|A|$ is the number of actions. The $\epsilon$-greedy policy is an example of an exploration-exploitation strategy, where the agent explores the environment by taking random actions (exploration) while also exploiting the current knowledge of the environment by taking the best-known action (exploitation). The value of $\epsilon$ is typically decayed over time, so that the agent explores more in the beginning and exploits more later on.

Monte Carlo methods are model-free, meaning they do not require a model of the environment (transition probabilities and rewards). However, they can be inefficient due to high variance and the need for many episodes to achieve accurate value estimates. They also require careful handling of the exploration-exploitation trade-off. The two main approaches for Monte Carlo methods are first-visit and every-visit methods.

First-Visit MC: In this approach, the return for a state is averaged over all first visits to that state in each episode.
Every-Visit Monte Carlo: Here, the return is averaged over every visit to the state, not just the first visit in each episode.

Monte Carlo Policy Iteration involves alternating between policy evaluation (estimating the value function of the current policy using Monte Carlo methods) and policy improvement (improving the policy based on the current value function estimate). This process is repeated until the policy converges to the optimal policy. Policy iteration with exact evaluation is guaranteed to do so in a finite number of steps for finite MDPs (Puterman 2014); with Monte Carlo estimates the evaluation is noisy and this guarantee does not carry over.

To find the optimal policy, a balance between exploration and exploitation must be maintained. This is achieved through strategies like $\epsilon$-greedy exploration. In Monte Carlo Control, the policy is often improved in a greedy manner based on the current value function estimate.

Recall that an arithmetic average can be updated recursively \[ \bar{x}_n = \frac{1}{n}\sum_{i=1}^n x_i = \frac{1}{n}\left(x_n + \sum_{i=1}^{n-1} x_i\right) = \frac{1}{n}\left(x_n + (n-1)\bar{x}_{n-1}\right) = \bar{x}_{n-1} + \frac{1}{n}(x_n - \bar{x}_{n-1}). \] This is called a running average. We can use this recursion to update the value function $V(s)$ incrementally, each time we visit state $s$ at time $t$. \[ V(s_t) = V(s_t) + \frac{1}{N(s_t)}(G_t - V(s_t)), \] where $N(s_t)$ is the number of times we visited state $s_t$ up to and including time $t$ and $G_t$ is the return at time $t$. This is called every-visit Monte Carlo method. Alternatively, we can use first-visit Monte Carlo method, where we update the value function only at the first visit to state $s$ in each episode, with $N(s)$ counting only those first visits.

Alternatively, we can use a learning rate $\alpha$ \[ V_{n+1} = V_n + \alpha(G_n - V_n). \] This is called a constant step size update. The learning rate is a hyperparameter that needs to be tuned. The constant step size update is more convenient because it does not require keeping track of the number of visits to each state. The constant step size update is also more robust to non-stationary problems.

Temporal Difference Learning (TD Learning) Similar to MC, TD methods learn directly from raw experience without a model of the environment. However, unlike MC methods, TD methods update value estimates based on other learned estimates, without waiting for the end of an episode. This is called bootstrapping. TD methods combine the sampling of Monte Carlo methods with the bootstrapping of dynamic programming. However, they are more complex than MC methods and require careful tuning of the learning rate.

A simple TD method is TD(0), which updates value estimates based on the current reward and the estimated value of the next state. The update rule is \[ V(S_t) = V(S_t) + \alpha(R_{t+1} + \gamma V(S_{t+1}) - V(S_t)), \] where $\alpha$ is the learning rate. The TD(0) method is also called one-step TD because it only looks one step ahead. The $R_{t+1} + \gamma V(S_{t+1})$ term is called the TD target and is a biased estimate of $V(S_t)$. The difference $R_{t+1} + \gamma V(S_{t+1}) - V(S_t)$ is called the TD error. The TD target is an estimate of the return $G_t$ and the TD error is the difference between the TD target and the current estimate $V(S_t)$. Although TD algorithms have lower variance than MC methods, they have higher bias. In practice, TD methods have usually been found to converge faster than constant-step Monte Carlo on stochastic tasks (Sutton and Barto 2018).

9.4 Q-Learning

Q-learning is an off-policy algorithm that learns the optimal policy by directly estimating the optimal action-value function $Q^*(s,a)$. The algorithm iteratively updates the action-value function using the Bellman optimality backup. Off-policy means that the algorithm learns the optimal policy while following a different policy. The algorithm can learn the optimal policy while following a random policy, for example, provided that every state-action pair is visited infinitely often and the learning rates satisfy $\sum_t \alpha_t = \infty$ and $\sum_t \alpha_t^2 < \infty$ (Watkins and Dayan 1992). The update rule is \[ Q(S_t,A_t) = Q(S_t,A_t) + \alpha(R_{t+1} + \gamma \max_a Q(S_{t+1},a) - Q(S_t,A_t)), \] where $\alpha$ is the learning rate. We can simplify the update rule to \[ Q(S_t,A_t) = (1-\alpha)Q(S_t,A_t) + \alpha(R_{t+1} + \gamma \max_a Q(S_{t+1},a)). \] At a terminal state $S_{t+1}$ the maximum is taken to be zero.

Initialize $Q(s,a)$ arbitrarily
Repeat for each episode:
1. Initialize $S$
2. Repeat for each step of the episode:
  1. Choose $A$ from $S$ using policy derived from $Q$ (e.g., $\epsilon$-greedy)
  2. Take action $A$, observe $R$, $S'$
  3. $Q(S,A) = Q(S,A) + \alpha(R + \gamma \max_a Q(S',a) - Q(S,A))$
  4. $S = S'$
3. Until $S$ is terminal

9.5 Bayesian Optimization

Bayesian Optimization solves an optimization problem using sequential design of experiments, which is the same in spirit as reinforcement learning but focuses on finding the optimum of a static function rather than a policy.

Bayesian optimization is a sequential design strategy for global optimization of black-box functions that does not assume any functional forms. It is particularly useful when the objective function is expensive to evaluate. Bayesian optimization uses a surrogate model to approximate the objective function and an acquisition function to decide where to sample next. The surrogate model is typically a Gaussian process (GP) model, which is a probabilistic model that defines a distribution over functions. The acquisition function is a heuristic that trades off exploration and exploitation to decide where to sample next. Bayesian optimization is derivative-free: it needs only evaluations of $f$, never its gradient. It targets the global optimum rather than a local one. Unlike convex optimization, it offers no guarantee of finding the true global optimum; instead, it seeks near-global optima efficiently by balancing exploration and exploitation. It is sample-efficient, often finding good solutions with far fewer function evaluations than grid search or random search (Jones, Schonlau, and Welch 1998). However, it can be slow in practice and is not suitable for high-dimensional problems.

Given a function $f(x)$ that is not known analytically (it can represent, for example, the output of a complex computer program), the goal is to optimize \[ x^* = \arg\min_x f(x). \]

The Bayesian approach to this problem is the following:

Define a prior distribution over $f(x)$
Calculate $f$ at a few points $x_1, \ldots, x_n$
Repeat until convergence:
1. Update the prior to get the posterior distribution over $f(x)$
2. Choose the next point $x^+$ to evaluate $f(x)$
3. Calculate $f(x^+)$
Pick $x^*$ that corresponds to the smallest value of $f(x)$ among evaluated points

The prior distribution is typically a Gaussian process (GP). The GP model is defined by a mean function $m(x)$ and a covariance function $k(x,x')$. The mean function is typically set to zero. The covariance function is typically a squared exponential function \[ k(x,x') = \sigma_f^2 \exp\left(-\frac{(x-x')^2}{2l^2}\right), \] where $\sigma_f^2$ is the signal variance, and $l$ is the length scale. The covariance function defines the similarity between two points $x$ and $x'$. The covariance function is also called a kernel function.

Now we need to decide where to sample next, which is the job of the acquisition function. The expected improvement (EI) function is a popular acquisition function (Jones, Schonlau, and Welch 1998). Suppose \[ f^* = \min_i y_i \] is the smallest response among the evaluated points $y_1,\ldots,y_n$. At a candidate point $x$, let $Y(x)$ denote the posterior predictive for the response there, which is the random quantity we have not yet observed. The expected improvement function is defined as \[ a(x) = \E{\max(0, f^* - Y(x))}. \] The function that we calculate the expectation of \[ u(x) = \max(0, f^* - Y(x)) \] is the utility function. Thus, the acquisition function is the expected value of the utility function.

The acquisition function is high when $Y(x)$ is likely to be lower than $f^*$, and low when it is likely to be higher. Given the GP prior, we can calculate the acquisition function analytically. The posterior predictive is normal, $Y(x) \sim N(\mu,\sigma^2)$ with $\mu = \mu(x)$ and $\sigma^2 = \sigma^2(x)$, so writing $y$ for the dummy integration variable the acquisition function is \[\begin{align*} a(x) &= \E{\max(0, f^* - Y(x))} \\ &= \int_{-\infty}^{\infty} \max(0, f^* - y) \phi(y,\mu,\sigma^2) dy \\ & \text{Since we are interested in improvement } f^* - y > 0 \text{, i.e., } y < f^* \text{,} \\ & \text{the integral is from } -\infty \text{ to } f^* \\ &= \int_{-\infty}^{f^*} (f^* - y) \phi(y,\mu,\sigma^2) dy \end{align*}\] where $\phi(y,\mu,\sigma^2)$ is the probability density function of the normal distribution. A useful identity is \[ \int y \phi(y,\mu,\sigma^2) dy =\frac{1}{2} \mu ~ \text{erf}\left(\frac{y-\mu }{\sqrt{2} \sigma }\right)-\frac{\sigma e^{-\frac{(y-\mu )^2}{2 \sigma ^2}}}{\sqrt{2 \pi }}, \] where erf is the error function. Thus, \[ \int_{-\infty}^{f^*} y \phi(y,\mu,\sigma^2) dy = \frac{1}{2} \mu (1+\text{erf}\left(\frac{f^*-\mu }{\sqrt{2} \sigma }\right))-\frac{\sigma e^{-\frac{(f^*-\mu )^2}{2 \sigma ^2}}}{\sqrt{2 \pi}} = \mu \Phi(f^*,\mu,\sigma^2) - \sigma^2 \phi(f^*,\mu,\sigma^2), \] where $\Phi(y,\mu,\sigma^2)$ is the cumulative distribution function of the normal distribution. We can write the acquisition function as \[ a(x) = \sigma^2 \phi(f^*,\mu,\sigma^2) + (f^*-\mu)\Phi(f^*,\mu,\sigma^2) \]

We can implement it

# p must be the predGP output on the grid the caller will index into.
acq <- function(p, fstar) {
  d <- fstar - p$mean
  s <- sqrt(diag(p$Sigma))
  return(s * dnorm(d / s) + d * pnorm(d / s))
}

Example 9.7 (Taxi Fleet Optimisation) We will use the taxi fleet simulator from Emukit project (Paleyes, Mahsereci, and Lawrence 2023). For a given demand (the frequency of trip requests) and the number of taxis in the fleet, it simulates the taxi fleet operations and calculates the profit. The simulator is a black-box function, meaning it does not have an analytical form and can only be evaluated at specific points. The goal is to find the optimal number of taxis in the fleet that maximizes the profit. We will use Bayesian optimization to solve this problem.

We start with an initial set of three designs $x = (10,90,30)$, where $x$ is the number of taxis in the fleet, and observe the corresponding profits $(3.1, 3.6, 6.6)$. When $x=10$, the demand for taxis exceeds the supply and passengers need to wait for their rides, leading to missed profit opportunities. At another extreme, when we have 90 taxis, the profit is slightly better. However, there are many empty taxis, which is not profitable. The optimal number of taxis must be somewhere in the middle. Finally, we try 30 taxis and observe that the profit is higher than both of our previous attempts. However, should we increase or decrease the number of taxis from here? We can use Bayesian optimization to answer this question. First, we use a convenience function, plotgp, that plots the GP emulator: the observations, the posterior median and a 90% band.

Now, we fit the GP emulator using our initial set of observed taxi-profit pairs, the result is shown in Figure 9.7.

library(laGP)
x <- matrix(c(10, 90, 30), ncol = 1)
xx <- seq(1, 100, length = 500)
XX <- matrix(xx, ncol = ncol(x))
profit <- -c(3.1, 3.6, 6.6)
gp <- newGP(x, profit, 1000, 1e-6, dK = TRUE)
p <- predGP(gp, XX)
plotgp(x, profit, XX, p)

Figure 9.7: Initial GP emulator fit to the taxi-profit pairs.

Instead of maximizing the profit, we minimize the negative profit. We see that there is potentially a better value at around 50 taxis. We can use the acquisition function to decide where to sample next. We define two functions: nextsample that uses the acquisition function to decide where to sample next and updgp that updates the GP emulator with the new sample. Then we call those two functions four times. The acquisition function successively suggests 48, 58, 45, and 100 taxis. The fleet simulator was run at 44, 57, 45, and 100, so the first two evaluations sit a few taxis below the acquisition maximum; this shifts the design slightly but not the region the search settles on.

nextsample <- function() {
  ei <- acq(p, min(profit))
  plot(xx, ei, type = "l", col = 2, lwd = 2, xlab = "Number of taxis", ylab = "Expected improvement")
  xnext <- as.integer(xx[which.max(ei)])
  return(xnext)
}
updgp <- function(xnext, f) {
  profit <<- c(profit, f)
  x <<- rbind(x, xnext)
  gp <<- newGP(x, profit, 1000, 1e-6, dK = TRUE)
  p <<- predGP(gp, XX)
  plotgp(x, profit, XX, p) + theme_gray(base_size = 8)
}
invisible(nextsample()) # suggests 48; simulator was run at 44
updgp(44, -8.4)
invisible(nextsample()) # suggests 58; simulator was run at 57
updgp(57, -7.1)
invisible(nextsample()) # suggests 45
updgp(45, -8.5)
invisible(nextsample()) # suggests 100
updgp(100, -3.3)

If we run nextsample one more time, we get 47, close to our current best of 45. Further, the model is confident at this location. It suggests that further evaluations would add little, and we stop here, as shown in Figure 9.9.

nextsample()
## 47

Figure 9.9: The acquisition function after four optimization steps, with seven design points. The maximum is at 47.

9.6 Concluding Remarks

This chapter developed three interconnected frameworks for sequential decision-making under uncertainty: multi-armed bandits for experimentation, MDPs and Q-learning for dynamic optimization, and Bayesian optimization for expensive black-box functions. The common thread is the exploration-exploitation tradeoff, balancing information gathering against immediate reward. Thompson sampling and expected improvement each resolve this tradeoff by maintaining and updating a probabilistic model of the unknown. Bellman’s principle of optimality reduces a sequential problem to one-step problems, and Q-learning applies it without a model, from sampled experience.

Exercises

Pen-and-Pencil

Exercise 9.1 (Bellman Equation) Consider a simple grid world where a robot can move in 4 directions. The state is the robot’s position, and the reward is $-1$ for each step, including the step into the goal (to encourage efficiency), and the goal is terminal with value 0.

In a 3×3 grid:

+---+---+---+
| 1 | 2 | 3 |
+---+---+---+
| 4 | 5 | 6 |
+---+---+---+
| 7 | 8 | G |
+---+---+---+

State 9 is the goal (G). Discount factor $\gamma = 0.9$.

Write the Bellman optimality equation for this problem.
What is $V^*(G) = V^*(9)$?
Assuming deterministic transitions, calculate $V^*(8)$ and $V^*(6)$.
Calculate $V^*(5)$.

Exercise 9.2 (Secretary Problem) You are hiring for a position and will interview exactly 10 candidates. After each interview, you must immediately decide to hire or reject (no callbacks). Candidates can be ranked relatively but not absolutely. You want to maximize the probability of hiring the best candidate.

What is the optimal strategy? (Hint: reject the first $r$ candidates, then hire the first one better than all previous)
For $n = 10$ candidates, what is the optimal $r$?
What is the probability of success with the optimal strategy?
(Optional, with a computer.) Simulate this problem 10,000 times to empirically verify your answer.

Exercise 9.3 (Q-Learning Update) A simple MDP has states {A, B, C} and actions {left, right}. The discount factor is $\gamma = 0.9$ and learning rate is $\alpha = 0.1$.

Current Q-table:

	left	right
A	2.0	3.0
B	1.5	2.5
C	0.0	0.0

The agent is in state A, takes action “right”, receives reward $r = 5$, and transitions to state B.

Write the Q-learning update formula.
Calculate the new value of $Q(A, \text{right})$.
If instead the agent transitioned to state C (terminal), what would be the update?
Suppose the environment is deterministic: A with action right goes to B with reward 5, B with action right goes to C with reward 2, B with action left goes to C with reward 1, and C is terminal. To what value does Q(A, right) converge under Q-learning, if every state-action pair keeps being visited?

Exercise 9.4 (TD(0) versus Monte Carlo) Consider a 3-state Markov reward chain $1 \to 2 \to 3$, where state 3 is terminal. The reward for transitioning $1\to 2$ is $R_1$ and for $2 \to 3$ is $R_2$, both random with means $\E{R_1}=3$ and $\E{R_2}=5$. Discount factor $\gamma = 0.9$, learning rate $\alpha = 0.1$. Current value estimates are $V(1)=2.0$, $V(2)=4.0$, $V(3)=0$ (terminal).

The agent is in state 1, observes reward $r=3$, and transitions to state 2. Compute the TD(0) update \[ V(1) \leftarrow V(1) + \alpha\big[r + \gamma V(2) - V(1)\big]. \]
In the same episode, the agent then observes reward $r_2=5$ transitioning from state 2 into the terminal state 3. Using the full realized return $G_1 = r + \gamma r_2$ from state 1, compute the every-visit Monte Carlo update $V(1) \leftarrow V(1) + \alpha(G_1-V(1))$, and compare to (a).
(Optional, with a computer.) Simulate 2,000 independent replications of learning $V(1)$ from only 10 episodes each, with $R_1 \sim N(3,3^2)$ and $R_2 \sim N(5,3^2)$ resampled every episode, and both estimators starting from $V(1)=V(2)=0$. Compare the bias (distance of the mean estimate from the true $V(1)$) and the variance of the TD(0) and Monte Carlo estimators of $V(1)$ after 10 episodes.

Exercise 9.5 (Regret in Bandits) A 3-armed bandit has true success probabilities:

Arm 1: $\theta_1 = 0.3$
Arm 2: $\theta_2 = 0.5$ (optimal)
Arm 3: $\theta_3 = 0.4$

An agent uses $\epsilon$-greedy with $\epsilon = 0.1$, playing 1000 rounds.

What is the expected regret per round of pulling arm 1 instead of arm 2?
What is the expected regret per round of the $\epsilon$-greedy strategy (assuming it correctly identifies arm 2 as best)?
Calculate the expected cumulative regret after 1000 rounds.
(Optional, with a computer.) Compare to Thompson Sampling: simulate both algorithms and plot cumulative regret.

Exercise 9.6 (True/False: Reinforcement Learning)

The Bellman equation is always linear in the value function.
Q-learning is an off-policy algorithm.
In the secretary problem, the probability of selecting the best candidate approaches 1 as the number of candidates increases.
Thompson Sampling requires knowing the true reward distributions of each arm.
In a continuing infinite-horizon MDP in which rewards stay bounded away from zero, the discount factor $\gamma$ must be strictly less than 1 for the value function to be finite.

Computing

Exercise 9.7 (Multi-Armed Bandit: Thompson Sampling) A website tests 3 button colors: Red (A), Green (B), Blue (C). Each has an unknown click-through rate $\theta_i$. After some initial testing, the posterior distributions are:

Red: Beta(15, 85)
Green: Beta(22, 78)
Blue: Beta(18, 82)

Calculate the posterior mean click-through rate for each button.
Using Thompson Sampling, simulate which button would be selected in the next trial.
Repeat part (b) 1000 times and report the selection frequencies.
After 100 more trials with 25 clicks on Green, update the Green posterior and recalculate.

Exercise 9.8 (Multi-Armed Bandit A/B Test) Consider a complex A/B experiment, with 6 alternatives, when you have 5 variations to your page, plus the original.

Use Bonferroni correction for multiple comparisons. Calculate the significance level for each of the 5 tests and find the number of samples per arm needed to achieve a power of 0.95, in a two-sided test, against an improvement from the original’s 4% conversion rate to 5%. Assume that the overall significance level is .05.
Implement TS for the same experiment, with 100 visits per day. Assume an original arm with a 4% conversion rate, and an optimal arm with a 5% conversion rate. The other 4 arms include one suboptimal arm that beats the original with a conversion rate of 4.5%, and three inferior arms with rates of 3%, 2%, and 3.5%. End the experiment when the potential value remaining, the stopping rule of this chapter, falls below 1% of the champion’s conversion rate, or after 919 days, the length of the classical test, whichever comes first. Plot the savings from a six-armed experiment, relative to a Bonferroni adjusted power calculation for a classical experiment. First plot should show the number of days required to end the experiment, with the vertical line showing the time required by the classical power calculation. The second plot should show the number of conversions that were saved by the bandit. What is the overall cost savings due to ending the experiment more quickly, and due to the experiment being less wasteful while it is running?
Run your simulator 500 times and show the history of the serving weights for all the arms in the first of our 500 simulation runs. Comment on the results.
Plot the daily cost of running the multi-armed bandit relative to an “oracle” strategy of always playing arm 2, the optimal arm.

Exercise 9.9 (Contextual Bandits: Weekday versus Weekend) The chapter’s Contextual Bandits section describes, without ever simulating it, an experiment in which “Arm A performs slightly better during the weekdays when users browse but don’t make purchases, while Arm B excels during the weekends when users are more likely to make purchases.” Make this concrete. Let $x_t=1$ denote a weekend, with $\prob{x_t=1}=2/7$, and $x_t=0$ a weekday. The true click-through rates are: arm A is $0.06$ on weekdays and $0.03$ on weekends; arm B is $0.03$ on weekdays and $0.08$ on weekends. A fresh context $x_t$ is drawn each round.

Implement a context-blind Thompson sampling bandit that keeps a single $Beta(1,1)$ posterior per arm (as in the chapter’s bandit() function) and never observes $x_t$.
Implement a context-aware Thompson sampling bandit that keeps a separate $Beta(1,1)$ posterior for each of the four (arm, context) cells. Explain why this is the Bayesian-conjugate analogue of the chapter’s logistic bandit model $\mathrm{logit}\,\prob{\text{click}_a \mid \theta, x} = \beta_{0a} + \beta_a^\top x$ (an arm-by-context interaction), and why a version with one $\beta$ shared by the arms cannot represent a crossover in which different arms are optimal in different contexts.
Simulate $T=5000$ rounds for each policy and plot cumulative regret for both on the same axes, where the regret at round $t$ is $\max(\theta_A(x_t),\theta_B(x_t)) - \theta_{\text{chosen}}(x_t)$.
Report the final cumulative regret of each policy and their ratio. In terms of the blended (marginal) rates $\bar\theta_a = \prob{x=0}\theta_{a,\text{weekday}} + \prob{x=1}\theta_{a,\text{weekend}}$, explain which arm the context-blind bandit converges to playing almost exclusively, and why. By extending the simulation to $T=50{,}000$, argue whether the context-blind bandit’s regret grows linearly or sub-linearly in $T$, and contrast with the context-aware bandit and the logarithmic regret of Thompson sampling in the chapter.
Using only the (arm, context, outcome) triples logged by the context-blind bandit, fit glm(y ~ arm * ctx, family = binomial). Show that it recovers the true crossover even though the policy that generated the data never used context to choose an arm. Explain why this does not help the context-blind bandit’s own regret, and what it implies about needing context in the decision rule itself rather than only in a post-hoc analysis.

Exercise 9.10 (Bayesian Optimization: Expected Improvement) Continue the chapter’s Taxi Fleet Optimisation example (Example 9.7): a GP surrogate is fit to the initial design $x=(10,30,90)$ with observed profits $(3.1,6.6,3.6)$, using the chapter’s squared-exponential kernel and hyperparameters, on $y=-\text{profit}$ so that $x^*=\arg\min_x f(x)$.

Fit the GP to these three points and extract the posterior mean $\mu(x)$ and standard deviation $\sigma(x)$ at two candidate fleet sizes, $x=50$ and $x=70$.
By hand, i.e., directly from the closed-form formula rather than a library call, compute the Expected Improvement at each candidate, \[ z = \frac{f^*-\mu(x)}{\sigma(x)}, \qquad EI(x) = \big(f^*-\mu(x)\big)\Phi(z) + \sigma(x)\phi(z), \] where $f^*=\min(y_1,\ldots,y_n)$ is the best (most negative) value observed so far and $\phi,\Phi$ are the standard normal density and CDF. Which candidate does Bayesian optimization sample next, and why?
Since the actual taxi simulator is not available outside the Emukit demo, implement the full Bayesian optimization loop (fit GP, maximize $EI$ over a grid, evaluate, update) on the synthetic profit curve $\text{profit}(x)=9\exp\big(-(x-45)^2/1250\big)$, which peaks at essentially the same fleet size the chapter’s own search converged to. Starting from the same three points, run 5 iterations. Does the search find the true optimum $x^*=45$?

Agrawal, Shipra, and Navin Goyal. 2012. “Analysis of Thompson Sampling for the Multi-Armed Bandit Problem.” In Proceedings of the 25th Annual Conference on Learning Theory, 23:39.1–26. Proceedings of Machine Learning Research.

Bellman, Richard. 1957. Dynamic Programming. Princeton University Press.

Ferguson, Thomas S. 1989. “Who Solved the Secretary Problem?” Statistical Science 4 (3): 282–89. https://doi.org/10.1214/ss/1177012493.

Gilbert, John P., and Frederick Mosteller. 1966. “Recognizing the Maximum of a Sequence.” Journal of the American Statistical Association 61 (313): 35–73. https://doi.org/10.1080/01621459.1966.10502008.

Howard, Ronald A. 1960. Dynamic Programming and Markov Processes. MIT Press.

Jones, Donald R., Matthias Schonlau, and William J. Welch. 1998. “Efficient Global Optimization of Expensive Black-Box Functions.” Journal of Global Optimization 13 (4): 455–92. https://doi.org/10.1023/A:1008306431147.

Kaufmann, Emilie, Nathaniel Korda, and Rémi Munos. 2012. “Thompson Sampling: An Asymptotically Optimal Finite-Time Analysis.” In Algorithmic Learning Theory, 199–213. https://doi.org/10.1007/978-3-642-34106-9_18.

Paleyes, Andrei, Maren Mahsereci, and Neil D. Lawrence. 2023. “Emukit: A Python Toolkit for Decision Making Under Uncertainty.” In Proceedings of the 22nd Python in Science Conference, 68–75. https://doi.org/10.25080/gerudo-f2bc6f59-009.

Post, Thierry, Martijn J. van den Assem, Guido Baltussen, and Richard H. Thaler. 2008. “Deal or No Deal? Decision Making Under Risk in a Large-Payoff Game Show.” American Economic Review 98 (1): 38–71. https://doi.org/10.1257/aer.98.1.38.

Puterman, Martin L. 2014. Markov Decision Processes: Discrete Stochastic Dynamic Programming. John Wiley & Sons.

Scott, Steven L. 2015. “Multi-Armed Bandit Experiments in the Online Service Economy.” Applied Stochastic Models in Business and Industry 31 (1): 37–45.

Scott, Steven L. 2013. “Multi-Armed Bandit Experiments.” Google Analytics Content Experiments documentation. https://web.archive.org/web/20130704074743/https://support.google.com/analytics/answer/2844870.

Shannon, Claude E. 1951. “Presentation of a Maze-Solving Machine.” In Cybernetics: Transactions of the Eighth Conference. Josiah Macy Jr. Foundation.

Silver, David. 2015. “Lectures on Reinforcement Learning.” UCL course on reinforcement learning, Lecture 2: Markov Decision Processes. https://www.davidsilver.uk/teaching/.

Silver, David, Thomas Hubert, Julian Schrittwieser, Ioannis Antonoglou, Matthew Lai, Arthur Guez, Marc Lanctot, et al. 2017. “Mastering Chess and Shogi by Self-Play with a General Reinforcement Learning Algorithm.” arXiv. https://arxiv.org/abs/1712.01815.

Sutton, Richard S., and Andrew G. Barto. 2018. Reinforcement Learning: An Introduction. 2nd ed. MIT Press.

Thompson, William R. 1933. “On the Likelihood That One Unknown Probability Exceeds Another in View of the Evidence of Two Samples.” Biometrika 25 (3/4): 285–94.

Watkins, Christopher J. C. H., and Peter Dayan. 1992. “Q-Learning.” Machine Learning 8 (3–4): 279–92. https://doi.org/10.1007/BF00992698.