Markov Decision Processes and Dynamic Programming

An action can change both the immediate reward and the situations an agent will face next. A Markov decision process (MDP) describes that dependence. When its transitions and rewards are known, dynamic programming can compute a policy; when they are unknown, experience can supply information for learning. We will compare both approaches in the same small grid, then examine changes to the reward, transition model, and available data.

What the state must tell you

An MDP specifies states \(s\), available actions \(a\), transition probabilities \(P(s’\mid s,a)\), and expected immediate rewards \(r(s,a,s’)\). The Markov assumption says that the current state and action contain the information needed for the distribution of the next state and reward; additional history does not change that distribution. If location alone is insufficient because battery charge or an unseen opponent matters, the state representation or the model must change.

A policy is a decision rule. A deterministic policy \(\pi(s)\) chooses an action; a stochastic policy \(\pi(a\mid s)\) gives action probabilities. Thus a policy can itself be a function. Reinforcement learning need not receive correct-action labels, although demonstrations can also be used. Its reward signal evaluates consequences, which may include immediate and delayed effects.

For discount \(0\leq\gamma<1\), the return is \(G_t=\sum_{k=0}^{\infty}\gamma^k R_{t+k+1}\). The value \(V^\pi(s)=\mathbb E_\pi[G_t\mid S_t=s]\) describes following policy \(\pi\) from state \(s\). Action value \(Q^\pi(s,a)\) describes taking \(a\) first and following \(\pi\) afterward. In a finite MDP with bounded rewards, discounting makes these sums finite, even for policies that never terminate. A fixed-horizon or average-reward objective is a different formulation.

A Bellman backup combines one step with what follows

For a deterministic policy, the Bellman expectation equation is \[V^\pi(s)=\sum_{s’}P(s’\mid s,\pi(s))\left[r(s,\pi(s),s’)+\gamma V^\pi(s’)\right].\] It separates the next reward from the value of the remaining decisions. For a stochastic policy, average this expression over its action probabilities. The optimality equation replaces the fixed action with a maximum: \[V^*(s)=\max_a\sum_{s’}P(s’\mid s,a)\left[r(s,a,s’)+\gamma V^*(s’)\right].\] A backup evaluates this right-hand side using the current value estimates.

The grid has 25 states numbered row by row, starting at 0 in the upper-left corner. Actions 0–3 move up, down, left, and right; pushing against a boundary leaves the agent in place. Entering state 6 pays 5 and ends the episode; entering state 24 pays 40 and ends it. Entering states 11–13 costs 12; other nonterminal moves cost 1. Terminal states have zero subsequent reward and value. Rewards are paid on arrival, not repeatedly after termination.

For the two-step route to state 6, the return is \(-1+5\gamma\), or 1.5 at \(\gamma=0.5\). A penalty-free eight-step route to state 24 gives \(-\sum_{k=0}^{6}\gamma^k+40\gamma^7\). This route goes around the penalty strip. The printed undiscounted sum answers a different question from the discounted value, so both are labeled.

import numpy as np

H = W = 5
S = H * W
GOAL, NEAR = 24, 6
term = np.zeros(S, bool)
term[[GOAL, NEAR]] = True
active = np.flatnonzero(~term)
def step(s, a):
    if term[s]:
        return s
    row, col = divmod(s, W)
    dr, dc = [(-1,0), (1,0), (0,-1), (0,1)][a]
    return np.clip(row+dr, 0, H-1)*W + np.clip(col+dc, 0, W-1)
next_state = np.array([[step(s,a) for a in range(4)] for s in range(S)])
R = np.full(S, -1.0)
R[GOAL], R[NEAR] = 40.0, 5.0
R[[11,12,13]] = -12.0
reward = R[next_state]
reward[term] = 0.0

def action_values(V, gamma, rewards=reward):
    Q = np.zeros((S,4))
    Q[active] = rewards[active] + gamma * V[next_state[active]]
    return Q

def VI(gamma, rewards=reward, tol=1e-10):
    V = np.zeros(S)
    for sweep in range(1,10001):
        new = action_values(V,gamma,rewards).max(axis=1)
        if np.max(np.abs(new-V)) < tol:
            return new, sweep, sweep*len(active)*4
        V = new
    raise RuntimeError("Value iteration did not converge")

def evaluate(pi, gamma, rewards=reward):
    P = np.zeros((S,S))
    P[np.arange(S), next_state[np.arange(S),pi]] = 1.0
    return np.linalg.solve(np.eye(S)-gamma*P, rewards[np.arange(S),pi])

for gamma in (0.5,0.8,0.9,0.95,0.99):
    V,sweeps,_ = VI(gamma)
    pi = action_values(V,gamma).argmax(axis=1)
    state, total, discounted = 0, 0.0, 0.0
    path = [state]
    for t in range(60):
        if term[state]:
            break
        a = pi[state]
        total += reward[state,a]
        discounted += gamma**t * reward[state,a]
        state = next_state[state,a]
        path.append(int(state))
    print(f"gamma {gamma:.2f} sweeps {sweeps} V(start) {V[0]:.3f} "
          f"path return {discounted:.3f} undiscounted {total:.1f} end {state}")
    if gamma == 0.95:
        print("path",path)
# gamma 0.50 sweeps 6 V(start) 1.500 path return 1.500 undiscounted 4.0 end 6
# gamma 0.80 sweeps 9 V(start) 4.437 path return 4.437 undiscounted 33.0 end 24
# gamma 0.90 sweeps 9 V(start) 13.915 path return 13.915 undiscounted 33.0 end 24
# gamma 0.95 sweeps 9 V(start) 21.900 path return 21.900 undiscounted 33.0 end 24
# path [0, 5, 10, 15, 20, 21, 22, 23, 24]
# gamma 0.99 sweeps 9 V(start) 30.489 path return 30.489 undiscounted 33.0 end 24

At discount 0.5, the near terminal wins; at the other tested discounts, the far goal wins from state 0. These are solutions to different return objectives. The sum of discount weights is \(1/(1-\gamma)\), often called an effective horizon, but it is not a cutoff: rewards after that many steps still count. Discount can be selected for modeling or computational reasons; report the training discount and keep the evaluation objective fixed when comparing methods.

Value iteration repeatedly applies the optimality backup. For a finite discounted MDP, it contracts the maximum absolute value error by at most \(\gamma\) per sweep. Writing \(TV\) for an optimality backup, a small residual \(\|TV-V\|_\infty\) gives the error bound \(\|V-V^*\|_\infty\leq\|TV-V\|_\infty/(1-\gamma)\). A large discount can therefore make a given residual less reassuring. The code raises an error if the iteration limit is reached instead of silently declaring convergence.

Policy evaluation and policy improvement do different work

Policy iteration first evaluates its current policy, then improves it by choosing an action maximizing the one-step expression based on \(V^\pi\). With exact evaluation, that improvement cannot reduce value. For a fixed policy, the values solve \((I-\gamma P_\pi)V=r_\pi\). The helper evaluate solves this linear system; iterative policy evaluation instead repeatedly applies the expectation backup. The method below uses iterative evaluation restarted from zero at each outer iteration.

def PI(gamma, tol=1e-10):
    pi = np.zeros(S,int)
    terms = 0
    for outer in range(1,201):
        V = np.zeros(S)
        for _ in range(10000):
            new = np.zeros(S)
            new[active] = reward[active,pi[active]] + gamma*V[next_state[active,pi[active]]]
            terms += len(active)
            if np.max(np.abs(new-V)) < tol:
                V = new
                break
            V = new
        else:
            raise RuntimeError("Policy evaluation did not converge")
        Q = action_values(V,gamma)
        terms += 4*len(active)
        candidate = Q.argmax(axis=1)
        keep = Q[np.arange(S),pi] >= Q.max(axis=1)-1e-10
        candidate[keep] = pi[keep]
        candidate[term] = 0
        if np.array_equal(candidate,pi):
            return V,pi,outer,terms
        pi = candidate
    raise RuntimeError("Policy iteration did not stabilize")

print("gamma VI sweeps action terms PI rounds action terms max value gap")
for gamma in (0.8,0.95,0.99):
    v1,it1,n1 = VI(gamma)
    pi1 = action_values(v1,gamma).argmax(axis=1)
    n1 += 4*len(active)
    v2,pi2,it2,n2 = PI(gamma)
    exact = evaluate(pi2,gamma)
    print(f"{gamma:.2f} {it1:9d} {n1:12d} {it2:9d} {n2:12d} {np.max(np.abs(v1-exact)):.2e}")
# gamma VI sweeps action terms PI rounds action terms max value gap
# 0.80         9          920         7        10810 0.00e+00
# 0.95         9          920         7        42596 0.00e+00
# 0.99         9          920         7       212152 0.00e+00

An “action term” here is one deterministic reward-plus-discounted-successor calculation. Value iteration uses four per nonterminal state; policy evaluation uses one, and improvement uses four. Terminal states contribute none. The VI count includes the final greedy-policy extraction, as well as the value sweeps. This is an operation count for these implementations, not elapsed time or a universal definition of one Bellman backup. The linear solve used to verify the final policy is outside the count.

Policy iteration takes fewer outer rounds here but performs many inner evaluation sweeps, especially for initial policies that loop and discounts near 1. Restarting evaluation from zero also costs work. Warm starts, linear-system solvers, sparse representations, or modified policy iteration can change the comparison. Fewer outer iterations alone does not identify the faster algorithm.

Modified policy iteration limits the evaluation work before improving again. In a scheme that chooses a greedy policy from the current values and then performs one expectation backup under that policy, the update equals a value-iteration sweep. More evaluation sweeps approach full policy evaluation. Neither endpoint is uniformly preferable. When several actions tie, retaining the current maximizing action avoids pointless policy changes; comparing final values also avoids treating different optimal policies as disagreement.

Q-learning estimates a backup from experience

Tabular Q-learning updates the experienced state–action entry: \[Q(s,a)\leftarrow Q(s,a)+\alpha\left[r+\gamma\max_b Q(s’,b)-Q(s,a)\right].\] The bracket is the temporal-difference error: a one-step target minus the old estimate. For a terminal successor, the continuation term is zero. For example, with old value 2, reward −1, successor maximum 4, discount 0.95, and step size 0.5, the target is 2.8 and the updated value is 2.4.

The behavior policy below is \(\varepsilon\)-greedy: with probability \(\varepsilon\) it draws a random action, otherwise it takes the current maximizing action. The target still maximizes over actions, rather than averaging over that behavior policy; this is the off-policy aspect of Q-learning. Only the chosen entry is directly updated in this table. The simulator generates transitions; the learner does not consult the transition table to evaluate untried actions.

Training starts at a uniformly sampled state, including terminals, and stops each rollout after at most 80 transitions. This artificial access improves coverage compared with always starting at state 0. The cap is a data-collection truncation, so the final nonterminal transition still bootstraps. If 80 steps were the actual task horizon, remaining time would need to enter the state and the terminal condition would change.

gamma = 0.95
V_star,_,_ = VI(gamma)
Q_star = action_values(V_star,gamma)
checkpoints = (100,500,2000,10000)

def qlearn_history(eps,seed):
    rng = np.random.default_rng(seed)
    Q = np.zeros((S,4))
    history = {}
    for ep in range(1,max(checkpoints)+1):
        state = int(rng.integers(S))
        for _ in range(80):
            if term[state]:
                break
            a = int(rng.integers(4)) if rng.random()<eps else int(Q[state].argmax())
            nxt = next_state[state,a]
            target = reward[state,a] + gamma*(0.0 if term[nxt] else Q[nxt].max())
            Q[state,a] += 0.5*(target-Q[state,a])
            state = nxt
        if ep in checkpoints:
            pi = Q.argmax(axis=1)
            optimal = np.isclose(Q_star[active,pi[active]],V_star[active],atol=1e-8,rtol=0).mean()
            history[ep] = (optimal,evaluate(pi,gamma)[0])
    return history

runs = {eps:[qlearn_history(eps,seed) for seed in range(3)] for eps in (0.0,0.1,0.3)}
print("episodes eps optimal-action fraction greedy V(start)")
for ep in checkpoints:
    for eps in runs:
        fraction,value = np.mean([run[ep] for run in runs[eps]],axis=0)
        print(f"{ep:8d} {eps:.1f} {fraction:23.4f} {value:15.3f}")
print(f"optimal V(start) {V_star[0]:.3f}")
# episodes eps optimal-action fraction greedy V(start)
#      100 0.0                  0.5217           3.750
#      100 0.1                  0.5797           3.750
#      100 0.3                  0.5362           3.750
#      500 0.0                  0.5217           3.750
#      500 0.1                  0.6522           3.750
#      500 0.3                  0.8696          21.900
#     2000 0.0                  0.5217           3.750
#     2000 0.1                  0.9420          21.900
#     2000 0.3                  1.0000          21.900
#    10000 0.0                  0.5217           3.750
#    10000 0.1                  1.0000          21.900
#    10000 0.3                  1.0000          21.900
# optimal V(start) 21.900

The table averages three seeds. The fraction counts any optimal action at each nonterminal state, including ties; the second number evaluates the learned greedy policy from state 0 under the original discounted objective. They measure different aspects of performance. A policy can be good from one start while acting poorly in states it never reaches from there.

At 10,000 episodes, the zero-exploration runs choose an optimal action in about 52% of nonterminal states and have value 3.750 from state 0; both exploring settings reach value 21.900. Greedy behavior can stop sampling useful alternatives, but these finite runs cannot show that it will fail forever. Zero initialization can itself encourage trying another action after negative rewards, and random restarts expose multiple states. Positive exploration improves coverage here; the table does not measure its reward cost during training or prove that one exploration rate is generally best.

The usual convergence result for stochastic tabular Q-learning requires repeated visits to every relevant state–action pair and per-pair step sizes satisfying \(\sum_t\alpha_t=\infty\) and \(\sum_t\alpha_t^2<\infty\), together with discounted finite-state assumptions and bounded rewards. The fixed step size 0.5 used in this deterministic demonstration is not that general convergence guarantee. Function approximation and offline data introduce further issues.

Exercises

1. Repeated bonuses and potential-based shaping. Compare a reward for entering selected cells with a shaping term that preserves the discounted objective.

import numpy as np

H = W = 5
S = H * W
GOAL, NEAR = 24, 6
term = np.zeros(S, bool)
term[[GOAL, NEAR]] = True
active = np.flatnonzero(~term)
def step(s, a):
    if term[s]:
        return s
    row, col = divmod(s, W)
    dr, dc = [(-1,0), (1,0), (0,-1), (0,1)][a]
    return np.clip(row+dr, 0, H-1)*W + np.clip(col+dc, 0, W-1)
next_state = np.array([[step(s,a) for a in range(4)] for s in range(S)])
R = np.full(S, -1.0)
R[GOAL], R[NEAR] = 40.0, 5.0
R[[11,12,13]] = -12.0
reward = R[next_state]
reward[term] = 0.0

def action_values(V, gamma, rewards=reward):
    Q = np.zeros((S,4))
    Q[active] = rewards[active] + gamma * V[next_state[active]]
    return Q

def VI(gamma, rewards=reward, tol=1e-10):
    V = np.zeros(S)
    for sweep in range(1,10001):
        new = action_values(V,gamma,rewards).max(axis=1)
        if np.max(np.abs(new-V)) < tol:
            return new, sweep, sweep*len(active)*4
        V = new
    raise RuntimeError("Value iteration did not converge")

def evaluate(pi, gamma, rewards=reward):
    P = np.zeros((S,S))
    P[np.arange(S), next_state[np.arange(S),pi]] = 1.0
    return np.linalg.solve(np.eye(S)-gamma*P, rewards[np.arange(S),pi])

gamma = 0.95
V_base,_,_ = VI(gamma)
for bonus in (0.0,0.5,1.0,2.0,5.0):
    shaped = reward.copy()
    shaped[active] += bonus*np.isin(next_state[active],[1,2,3,4])
    V,_,_ = VI(gamma,shaped)
    pi = action_values(V,gamma,shaped).argmax(axis=1)
    state,total = 0,0.0
    path = [0]
    for _ in range(80):
        if term[state]:
            break
        a = pi[state]
        total += reward[state,a]
        state = next_state[state,a]
        path.append(int(state))
    print(f"bonus {bonus:.1f} original discounted {evaluate(pi,gamma)[0]:.3f} "
          f"80-step undiscounted {total:.1f} terminal {bool(term[state])}")
    if bonus == 5:
        print("first states",path[:10])
phi = -np.array([abs(s//W-4)+abs(s%W-4) for s in range(S)],dtype=float)
phi[term] = 0
potential_reward = reward + gamma*phi[next_state]-phi[:,None]
v_shaped,_,_ = VI(gamma,potential_reward)
pi_shaped = action_values(v_shaped,gamma,potential_reward).argmax(axis=1)
print(f"potential shaping: max original value loss {np.max(V_base-evaluate(pi_shaped,gamma)):.2e}")
# bonus 0.0 original discounted 21.900 80-step undiscounted 33.0 terminal True
# bonus 0.5 original discounted 21.900 80-step undiscounted 33.0 terminal True
# bonus 1.0 original discounted 21.900 80-step undiscounted 33.0 terminal True
# bonus 2.0 original discounted 21.900 80-step undiscounted 33.0 terminal True
# bonus 5.0 original discounted -20.000 80-step undiscounted -80.0 terminal False
# first states [0, 1, 1, 1, 1, 1, 1, 1, 1, 1]
# potential shaping: max original value loss 0.00e+00

2. Plan for stochastic actions and evaluate the same objective. With probability “slip”, execute one of the other three actions uniformly. Compare policies using expected discounted return in that same environment.

import numpy as np

H = W = 5
S = H * W
GOAL, NEAR = 24, 6
term = np.zeros(S, bool)
term[[GOAL, NEAR]] = True
active = np.flatnonzero(~term)
def step(s, a):
    if term[s]:
        return s
    row, col = divmod(s, W)
    dr, dc = [(-1,0), (1,0), (0,-1), (0,1)][a]
    return np.clip(row+dr, 0, H-1)*W + np.clip(col+dc, 0, W-1)
next_state = np.array([[step(s,a) for a in range(4)] for s in range(S)])
R = np.full(S, -1.0)
R[GOAL], R[NEAR] = 40.0, 5.0
R[[11,12,13]] = -12.0
reward = R[next_state]
reward[term] = 0.0

def action_values(V, gamma, rewards=reward):
    Q = np.zeros((S,4))
    Q[active] = rewards[active] + gamma * V[next_state[active]]
    return Q

def VI(gamma, rewards=reward, tol=1e-10):
    V = np.zeros(S)
    for sweep in range(1,10001):
        new = action_values(V,gamma,rewards).max(axis=1)
        if np.max(np.abs(new-V)) < tol:
            return new, sweep, sweep*len(active)*4
        V = new
    raise RuntimeError("Value iteration did not converge")

def evaluate(pi, gamma, rewards=reward):
    P = np.zeros((S,S))
    P[np.arange(S), next_state[np.arange(S),pi]] = 1.0
    return np.linalg.solve(np.eye(S)-gamma*P, rewards[np.arange(S),pi])

gamma = 0.95

def slip_model(slip):
    P = np.zeros((S,4,S))
    r = np.zeros((S,4))
    for state in range(S):
        for a in range(4):
            for executed in range(4):
                prob = 1-slip if executed==a else slip/3
                P[state,a,next_state[state,executed]] += prob
                r[state,a] += prob*reward[state,executed]
    return P,r

def solve_model(P,r):
    V = np.zeros(S)
    for _ in range(10000):
        Q = r+gamma*np.einsum('sak,k->sa',P,V)
        new = Q.max(axis=1)
        if np.max(np.abs(new-V))<1e-11:
            return Q.argmax(axis=1)
        V = new
    raise RuntimeError("Stochastic value iteration did not converge")

def policy_values(pi,P,r):
    Ppi = P[np.arange(S),pi]
    rpi = r[np.arange(S),pi]
    discounted = np.linalg.solve(np.eye(S)-gamma*Ppi,rpi)[0]
    finite = np.zeros(S)
    for _ in range(80):
        finite = rpi+Ppi@finite
    return discounted,finite[0]

pi0 = solve_model(*slip_model(0.0))
print("slip discounted aware unaware gain | 80-step undiscounted aware unaware")
for slip in (0.0,0.1,0.3,0.5):
    P,r = slip_model(slip)
    pi = solve_model(P,r)
    va,ua = policy_values(pi,P,r)
    vd,ud = policy_values(pi0,P,r)
    print(f"{slip:.1f} {va:16.4f} {vd:7.4f} {va-vd:7.4f} | {ua:10.4f} {ud:7.4f}")
# slip discounted aware unaware gain | 80-step undiscounted aware unaware
# 0.0          21.9002 21.9002  0.0000 |    33.0000 33.0000
# 0.1          18.6610 18.5591  0.1020 |    30.2880 30.2386
# 0.3          10.0823  9.5654  0.5169 |    21.6573 21.1807
# 0.5           1.2899 -1.8528  3.1427 |     3.2148  4.2956

3. Fit values from a fixed log. Hold unseen actions at different values, and compare unconstrained maximization with maximization restricted to logged actions.

import numpy as np

H = W = 5
S = H * W
GOAL, NEAR = 24, 6
term = np.zeros(S, bool)
term[[GOAL, NEAR]] = True
active = np.flatnonzero(~term)
def step(s, a):
    if term[s]:
        return s
    row, col = divmod(s, W)
    dr, dc = [(-1,0), (1,0), (0,-1), (0,1)][a]
    return np.clip(row+dr, 0, H-1)*W + np.clip(col+dc, 0, W-1)
next_state = np.array([[step(s,a) for a in range(4)] for s in range(S)])
R = np.full(S, -1.0)
R[GOAL], R[NEAR] = 40.0, 5.0
R[[11,12,13]] = -12.0
reward = R[next_state]
reward[term] = 0.0

def action_values(V, gamma, rewards=reward):
    Q = np.zeros((S,4))
    Q[active] = rewards[active] + gamma * V[next_state[active]]
    return Q

def VI(gamma, rewards=reward, tol=1e-10):
    V = np.zeros(S)
    for sweep in range(1,10001):
        new = action_values(V,gamma,rewards).max(axis=1)
        if np.max(np.abs(new-V)) < tol:
            return new, sweep, sweep*len(active)*4
        V = new
    raise RuntimeError("Value iteration did not converge")

def evaluate(pi, gamma, rewards=reward):
    P = np.zeros((S,S))
    P[np.arange(S), next_state[np.arange(S),pi]] = 1.0
    return np.linalg.solve(np.eye(S)-gamma*P, rewards[np.arange(S),pi])

gamma = 0.95
V_star,_,_ = VI(gamma)
Q_star = action_values(V_star,gamma)
pi_star = Q_star.argmax(axis=1)

def collect(eps,episodes=400):
    rng = np.random.default_rng(0)
    data = {}
    for _ in range(episodes):
        state = int(rng.integers(S))
        for _ in range(80):
            if term[state]:
                break
            a = int(rng.integers(4)) if rng.random()<eps else int(pi_star[state])
            nxt = next_state[state,a]
            data[state,a] = (reward[state,a],nxt)
            state = nxt
    return data

def offline(data,unseen_value,restrict=False):
    seen = np.zeros((S,4),bool)
    logged_r = np.zeros((S,4))
    logged_next = np.zeros((S,4),int)
    for (state,a),(r,nxt) in data.items():
        seen[state,a] = True
        logged_r[state,a],logged_next[state,a] = r,nxt
    if restrict and not np.all(seen[active].any(axis=1)):
        raise ValueError("This example requires an action logged in every nonterminal state")
    Q = np.full((S,4),float(unseen_value))
    Q[seen] = 0.0
    Q[term] = 0.0
    for _ in range(2000):
        choices = np.where(seen,Q,-np.inf) if restrict else Q
        values = choices.max(axis=1)
        values[term] = 0.0
        new = Q.copy()
        new[seen] = logged_r[seen]+gamma*values[logged_next[seen]]
        if np.max(np.abs(new-Q))<1e-10:
            Q = new
            break
        Q = new
    else:
        raise RuntimeError("Offline backups did not converge")
    if restrict:
        Q = np.where(seen,Q,-np.inf)
    return Q.argmax(axis=1)

print("behavior eps pairs coverage mode optimal-action fraction V(start)")
for eps in (0.0,0.05,0.2,0.5,1.0):
    data = collect(eps)
    for name,initial,restrict in (("unseen 0",0,False),("unseen 50",50,False),("logged only",0,True)):
        pi = offline(data,initial,restrict)
        fraction = np.isclose(Q_star[active,pi[active]],V_star[active],rtol=0,atol=1e-8).mean()
        print(f"{eps:12.2f} {len(data):5d} {len(data)/(4*len(active)):8.4f} "
              f"{name:11s} {fraction:23.4f} {evaluate(pi,gamma)[0]:8.3f}")
# behavior eps pairs coverage mode optimal-action fraction V(start)
#         0.00    23   0.2500 unseen 0                     1.0000   21.900
#         0.00    23   0.2500 unseen 50                    0.0000  -20.000
#         0.00    23   0.2500 logged only                  1.0000   21.900
#         0.05    60   0.6522 unseen 0                     1.0000   21.900
#         0.05    60   0.6522 unseen 50                    0.0435  -20.000
#         0.05    60   0.6522 logged only                  1.0000   21.900
#         0.20    86   0.9348 unseen 0                     1.0000   21.900
#         0.20    86   0.9348 unseen 50                    0.2174  -20.000
#         0.20    86   0.9348 logged only                  1.0000   21.900
#         0.50    92   1.0000 unseen 0                     1.0000   21.900
#         0.50    92   1.0000 unseen 50                    1.0000   21.900
#         0.50    92   1.0000 logged only                  1.0000   21.900
#         1.00    92   1.0000 unseen 0                     1.0000   21.900
#         1.00    92   1.0000 unseen 50                    1.0000   21.900
#         1.00    92   1.0000 logged only                  1.0000   21.900

Further reading

The next article, Deep Reinforcement Learning, extends these ideas to neural function approximation. For dynamic programming, see Berkeley’s policy-iteration notes. Ng, Harada, and Russell develop potential-based shaping, and Levine and colleagues’ offline RL tutorial explains learning from fixed datasets.