<a href="https://colab.research.google.com/github/probml/pyprobml/blob/master/notebooks/dtheory.ipynb" target="_parent"></a>

# Bayesian decision theory



```python
import numpy as np
np.set_printoptions(precision=3)

import pandas as pd
```

# From beliefs to actions

* Probabilistic inference let us compute our distribution over the possible **states of nature** $p(H=h|X=x)$; this is called our **belief state**. But ultimately we must convert beliefs into actions. 
* The optimal way to make **decisions under uncertainty** is to use **decision theory**.
* If the **decision maker** or **agent** picks action $a \in \cal{A}$ when the true state of nature is $h \in \cal{H}$, then we assume they incur a **loss** of $\ell(h,a)$.
* Since the state is usually hidden, we must compute the  the **expected loss** or **risk** for each action:
$$
R(a|x) = \sum_{h \in \cal{H}} p(h|x) \ell(h,a)
$$
* The **optimal policy** specifies which action to take for each possible observation so as to minimize this risk:
$$
a^* = \pi^*(x)  = \arg \min_{a \in \cal{A}} 
 \sum_{h \in \cal{H}} p(h|x) \ell(h,a)
$$
* Equivalently, we can **maximize expected utility**
$$
a^* = \pi^*(x) = \arg \max_{a \in \cal{A}} \sum_{h \in \cal{H}} p(h|x) U(h,a) 
$$


# Example:  Optimal policy for treating COVID patients

* Suppose we are a doctor who can either do nothing or give an expensive and painful drug to a patient.
* Suppose the patient either has or does not have COVID, and is either young or old.
* Suppose the loss function $\ell(h,a)$ is as shown below, where the units are [Quality Adjusted Life Years](https://en.wikipedia.org/wiki/Quality-adjusted_life_year#:~:text=The%20quality%2Dadjusted%20life%20year,one%20year%20in%20perfect%20health.). This matrix encodes the assumption that a young person is likely to live much longer than an old person, and that the drug is so bad that it is like losing 8 (quality-adjusted) years of life. (We will vary these assumptions later.)





* We assume the age is a visible variable, but the COVID status is hidden.
* The doctor can use Bayes rule to infer $p(h|x)$, where $x$ is a diagnostic test.
* Given the belief state, we can compute the optimal policy, which is shown below. 
* We see that we should only give the drug to young people who test positive.



## Sensitivity to assumptions

* If we reduce the cost of the drug from 8 to 5, we get the new optimal policy shown below.
* Now we see that we should give the drug to anyone who tests positive.



* If we keep the cost of the drug at 8, but increase the sensititivity of the test from 0.875 to 0.975, we get the new optimal policy shown below.
* Now we see that we should give the drug to anyone who tests positive, since the chance that they actually have COVID is higher.



## Code to reproduce the above.


```python
def make_loss_fun(loss_drug=15):
  loss_covid_young = 60
  loss_covid_old = 10
  loss_mat = np.zeros((4, 2))
  loss_mat[0, :] = [0, loss_drug]
  loss_mat[1, :] = [loss_covid_young, loss_drug]
  loss_mat[2, :] = [0, loss_drug]
  loss_mat[3, :] = [loss_covid_old, loss_drug]
  return loss_mat
```


```python
loss_mat = make_loss_fun()
records = []
records.append(('No COVID, young', loss_mat[0,0], loss_mat[0,1]))
records.append(('COVID, young', loss_mat[1,0], loss_mat[1,1]))
records.append(('No COVID, old', loss_mat[2,0], loss_mat[2,1]))
records.append(('COVID, old', loss_mat[3,0], loss_mat[3,1]))

df = pd.DataFrame.from_records(records, columns=['State', 'Nothing', 'Drugs'])
display(df)
```


<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>State</th>
      <th>Nothing</th>
      <th>Drugs</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>No COVID, young</td>
      <td>0.0</td>
      <td>15.0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>COVID, young</td>
      <td>60.0</td>
      <td>15.0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>No COVID, old</td>
      <td>0.0</td>
      <td>15.0</td>
    </tr>
    <tr>
      <th>3</th>
      <td>COVID, old</td>
      <td>10.0</td>
      <td>15.0</td>
    </tr>
  </tbody>
</table>
</div>



```python
def normalize(x):
  return x / np.sum(x)

def posterior_covid(observed, prevalence=None, sensitivity=None):
  # observed = 0 for negative test, 1 for positive test
  # hidden state = 0 if no-covid, 1 if have-covid
  if sensitivity is None:
    sensitivity = 0.875
  specificity = 0.975
  TPR = sensitivity; 
  FNR = 1-TPR
  TNR = specificity
  FPR = 1-TNR
  # likelihood(hidden, obs)
  likelihood_fn = np.array([[TNR, FPR], [FNR, TPR]])
  # prior(hidden)
  if prevalence is None:
    prevalence = 0.1
  prior = np.array([1-prevalence, prevalence])
  likelihood = likelihood_fn[:, observed].T
  posterior = normalize(prior * likelihood)
  return posterior
```


```python
def compute_belief(obs, prevalence=None, sensitivity=None):
  test_result, age = obs
  # returns bel[0:3] for states (not-covid+not-old, covid+not-old, not-covid+old, covid+old)
  post = posterior_covid(test_result, prevalence, sensitivity)
  bel = np.zeros(4)
  if age==0:
    bel[0:2] = post
  else:
    bel[2:4] = post
  return bel


def policy(obs):
  bel = compute_belief(obs)
  loss_mat = make_loss_fun(loss_drug)
  expected_loss = np.dot(bel, loss_mat)
  return np.argmin(expected_loss)


```


```python
def make_table(cost_drug, prevalence=None, sensitivity=None):
  records = []
  for test_result in [0,1]:
    for age in [0,1]:
      obs = (test_result, age)
      bel = compute_belief(obs, prevalence, sensitivity)    
      loss_mat = make_loss_fun(cost_drug)
      expected_loss = np.dot(bel, loss_mat)
      act = np.argmin(expected_loss)
      if age==0:
        post_covid = bel[1]
      else:
        post_covid = bel[3]
      records.append((test_result, age, post_covid, expected_loss[0], expected_loss[1], act))

  df = pd.DataFrame.from_records(records, columns=['test', 'age', 'pr(covid)', 'cost-noop', 'cost-drugs', 'action'])
  return df

```


```python
df = make_table(cost_drug=8)
display(df)
print(df.to_latex(index=False, float_format="%.2f" ))

```


<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>test</th>
      <th>age</th>
      <th>pr(covid)</th>
      <th>cost-noop</th>
      <th>cost-drugs</th>
      <th>action</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>0</td>
      <td>0.014045</td>
      <td>0.842697</td>
      <td>8.0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0</td>
      <td>1</td>
      <td>0.014045</td>
      <td>0.140449</td>
      <td>8.0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>1</td>
      <td>0</td>
      <td>0.795455</td>
      <td>47.727273</td>
      <td>8.0</td>
      <td>1</td>
    </tr>
    <tr>
      <th>3</th>
      <td>1</td>
      <td>1</td>
      <td>0.795455</td>
      <td>7.954545</td>
      <td>8.0</td>
      <td>0</td>
    </tr>
  </tbody>
</table>
</div>


    \begin{tabular}{rrrrrr}
    \toprule
     test &  age &  pr(covid) &  cost-noop &  cost-drugs &  action \\
    \midrule
        0 &    0 &       0.01 &       0.84 &        8.00 &       0 \\
        0 &    1 &       0.01 &       0.14 &        8.00 &       0 \\
        1 &    0 &       0.80 &      47.73 &        8.00 &       1 \\
        1 &    1 &       0.80 &       7.95 &        8.00 &       0 \\
    \bottomrule
    \end{tabular}
    



```python
df = make_table(cost_drug=5)
display(df)
print(df.to_latex(index=False, float_format="%.2f" ))
```


<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>test</th>
      <th>age</th>
      <th>pr(covid)</th>
      <th>cost-noop</th>
      <th>cost-drugs</th>
      <th>action</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>0</td>
      <td>0.014045</td>
      <td>0.842697</td>
      <td>5.0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0</td>
      <td>1</td>
      <td>0.014045</td>
      <td>0.140449</td>
      <td>5.0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>1</td>
      <td>0</td>
      <td>0.795455</td>
      <td>47.727273</td>
      <td>5.0</td>
      <td>1</td>
    </tr>
    <tr>
      <th>3</th>
      <td>1</td>
      <td>1</td>
      <td>0.795455</td>
      <td>7.954545</td>
      <td>5.0</td>
      <td>1</td>
    </tr>
  </tbody>
</table>
</div>


    \begin{tabular}{rrrrrr}
    \toprule
     test &  age &  pr(covid) &  cost-noop &  cost-drugs &  action \\
    \midrule
        0 &    0 &       0.01 &       0.84 &        5.00 &       0 \\
        0 &    1 &       0.01 &       0.14 &        5.00 &       0 \\
        1 &    0 &       0.80 &      47.73 &        5.00 &       1 \\
        1 &    1 &       0.80 &       7.95 &        5.00 &       1 \\
    \bottomrule
    \end{tabular}
    



```python
df = make_table(cost_drug=8, sensitivity=0.975)
display(df)
print(df.to_latex(index=False, float_format="%.2f" ))
```


<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>test</th>
      <th>age</th>
      <th>pr(covid)</th>
      <th>cost-noop</th>
      <th>cost-drugs</th>
      <th>action</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>0</td>
      <td>0.002841</td>
      <td>0.170455</td>
      <td>8.0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0</td>
      <td>1</td>
      <td>0.002841</td>
      <td>0.028409</td>
      <td>8.0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>1</td>
      <td>0</td>
      <td>0.812500</td>
      <td>48.750000</td>
      <td>8.0</td>
      <td>1</td>
    </tr>
    <tr>
      <th>3</th>
      <td>1</td>
      <td>1</td>
      <td>0.812500</td>
      <td>8.125000</td>
      <td>8.0</td>
      <td>1</td>
    </tr>
  </tbody>
</table>
</div>


    \begin{tabular}{rrrrrr}
    \toprule
     test &  age &  pr(covid) &  cost-noop &  cost-drugs &  action \\
    \midrule
        0 &    0 &       0.00 &       0.17 &        8.00 &       0 \\
        0 &    1 &       0.00 &       0.03 &        8.00 &       0 \\
        1 &    0 &       0.81 &      48.75 &        8.00 &       1 \\
        1 &    1 &       0.81 &       8.12 &        8.00 &       1 \\
    \bottomrule
    \end{tabular}
    


# Decision-theoretic classification

* In classification problems, the unknown state of nature is the "true" class label $y^* \in \cal{Y}$, where $\cal{Y} = \{1,\ldots,C\}$.
* We typically assume that the set of possible actions is to pick one of the labels, $\hat{y} \in \cal{Y}$. 
* We often assume a **zero-one** loss matrix $\ell(h,a)$, as follows:

$$
\begin{array}{c|cc}
 & \hat{y}=0 & \hat{y}=1  \\ \hline
 y^*=0 & 0 & 1  \\
 y^*=1 & 1 & 0
 \end{array}
$$
* We can write this more concisely as follows:
$$
\ell(y^*,\hat{y}) = I(y^* \neq \hat{y})
$$
* In this case,
the posterior expected loss is
$$
R(\hat{y}|x)
 = p(\hat{y} \neq y^* |x) = 1-p(y^*=\hat{y}|x)
$$
Hence the action that minimizes the expected loss is to choose the
most probable label:
$$
\pi(x)  = \arg \max_{y \in \cal{Y}} p(y|x)
$$
* This  corresponds to 
the **mode** of the posterior distribution,
also known as the **maximum a posteriori**
or  **MAP estimate**.



## Classification with a "reject" option

* In some applications (e.g., medicine, finance), we have an asymmetric loss, where the cost of a false positives is much higher than a cost of a false negative.
* In such settings, it is useful to allow for addition **rejection action**, in which we **refuse to classify** if we are uncertain. Such inputs can then be handled by some other system (e.g., a human).
* We define the loss function as follows:
$$
\ell(h, a) =
\left\{ \begin{array}{cc}
  0 & \mbox{if $h=a$ and $a \in \{1,\ldots,C\}$}\\
  \lambda_r & \mbox{if $a=0$}\\
  \lambda_e & \mbox{otherwise}
\end{array}
\right.
$$
where $\lambda_r$ is the cost of a rejection, and $\lambda_e$ is a cost of a misclassification error.
* It is easy to show that the optimal policy is as follows:
$$
\begin{align}
a^* &=
\begin{cases}
y^* & \text{ if } p(y^*|x) >  1-\frac{\lambda_r}{\lambda_e} \\
\text{reject} & \text{otherwise}
\end{cases}
\end{align}
$$
where $y^* = \arg \max_{y} p(y|x)$ is the MAP estimate.



# Bandit problems

## Multi-armed bandits

The classic multi-armed bandit is illustrated below.
The agent (the octopus) must choose which "arm" to pull.
Each slot machine (aka one-armed bandit) has a different, but unknown, payout rate (reward), denoted $R(a)$.





* In a MAB problem, the goal is to choose a sequence of actions (arm pulls) that maximizes the sum of expected rewards:
$$
a_{1:T} = \arg \max_{a_{1:T}} \sum_{t=}^T E_{a_t \sim \pi}[R(a_t | \theta)]
$$
where $T$ is the **finite horizon**, $\pi()$ is the (stochastic) **policy**, $R(a|\theta)$ is the reward function for action $a$ using model parameters $\theta$.
* In the **infinite horizon** case, we use a **discount factor** $\gamma < 1$ to ensure the sum converges:
$$
J = \sum_{t=1}^\infty \gamma^{t-1} E_{ a_t \sim \pi}[R(a_t | \theta)]
$$
* For binary rewards, we can use a **Bernoulli bandit**
$$
p(R|a) = \text{Ber}(R|\mu_a)
$$
where $\mu_a=E[R|a]=p(R=1|a)$.
* For real-valued rewards, we can use a **Gaussian bandit**
$$
p(R|a) = \text{Gauss}(R|\mu_a, \sigma_a^2)
$$


## Contextual bandits 

* In **contextual bandits**, the agent is in a random state $s_t$ at each step. States are drawn iid from some unknown distribution $p(s_t)$. The policy should now depend on the inputs, $a_t \sim \pi(s_t)$.
* The goal is to maximimize
$$
J = \sum_{t=1}^\infty \gamma^{t-1} E_{ s_t \sim p(s), a_t \sim \pi(a|s_t)}[R(s_t, a_t | \theta)]
$$
* There are numerous applications, e.g.
** **Online advertising**, where the action specifies which ad to show, and the state specifies user features (e.g., search history) and/or web page features (e.g., content).
** **Recommender systems**, where the action specifies which (set of) item(s)  to show, and the state specifies user features.
** **Clinical trials**, where the action specifies which drug to give, and the state specifies patient features.
* For binary rewards, we can use logistic regression model
$$
p(R|s,a;\theta) = \text{Ber}(R|\sigma(\theta^T \phi(s,a)))
$$
where $\phi(s,a)$ is a feature vector, possibly computed using a deep neural network, and $\sigma(x)=1/(1+e^{-x})$ is the sigmoid or logistic function.
* For real-valued rewards, we can replace the sigmoid output with a linear output, and use a Gaussian observation model.
* A key problem is that the rewards for each (state,action) pair (represented by model parameters $\theta$) are unknown.
* The agent can infer the parameters in an online (recursive) fashion given the data stream $D_{1:t} = (s_{1:t}, a_{1:t}, r_{1:t})$ using **sequential Bayesian updating**:
$$
\begin{align}
p(\theta|D_{1:t}) &= 
p(\theta|D_t, D_{1:t-1})
= \frac{p(D_t|\theta,D_{1:t-1}) p(\theta|D_{1:t-1})}
{p(D_t|D_{1:t-1})} 
\propto p(D_t|\theta) p(\theta|D_{1:t-1})
\end{align}
$$
* The likelihood is $p(D_t|\theta)=p(r_t|s_t, a_t,\theta)$ and the prior is the belief state from the previous time step, $p(\theta|D_{1:t-1})$.


## Exploration-exploitation tradeoff

* The agent must try each (state, action) pair (possibly multiple times) in order to collect data, but these experiments incur an opportunity cost. This is called the **exploration-exploitation tradeoff**.
* A common heuristic solution is the **upper confidence bound** method, which selects the action according to
$$
a_t = \arg \max_a E[R(s_t,a;\theta)] + c \;\text{std}[R(s_t,a;\theta)]
$$
* The  constant $c \geq 0$ reflects the amount of optimism about unknown rewards, and hence the amount of exploration.
* Below we illustrate the distribution over possible rewards for  a 3-armed bandit (without contextual features).



## Reinforcement learning


* In **reinforcement learning**, the setup is similar to a contextual bandit, except now the current state depends on the previously visited states, as well as the actions that the agent has taken. 
* If we assume a Markov model for the enviroment, we get a **Markov decision process** or **MDP**.
$$
J = \sum_{t=1}^\infty \gamma^{t-1} E_{ s_t \sim p(s|s_{t-1}, a_{t-1}; \theta), a_t \sim \pi(a|s_{1:t})}[R(s_t, a_t | \theta)]
$$
* Solving for the optimal policy when the parameters $\theta$ of the environment are unknown is challenging. This is beyond the scope of this tutorial.




```python

```

09951dfdd2b98b0a02e550e39ca8a3e20bb590f8
ipynb
Jupyter Notebook
notebooks/dtheory.ipynb
Prahitha/pyprobml
7dc7f58d4abd7a006e471bcaddd5dc24e294196f
MIT
notebooks/dtheory.ipynb
Prahitha/pyprobml
7dc7f58d4abd7a006e471bcaddd5dc24e294196f
MIT
2021-04-19T12:25:26.000Z
2021-04-19T12:25:26.000Z
notebooks/dtheory.ipynb
Nirzu97/pyprobml
397173b09668d21aac5c047b830db577cbb30530
MIT
Qwen/Qwen-72B
1. YES
2. YES