```python
%pylab inline
%load_ext autoreload
%autoreload 2
import matplotlib.pyplot as plt
import numpy as np
import numpy.random as npr
import sys
sys.path.append("../subjective-fits/")
import seaborn as sns
figsize(8,6)
```


```python
mpl.rcParams['xtick.labelsize'] = 22
mpl.rcParams['ytick.labelsize'] = 22
plt.rc('axes', labelsize=22)
plt.rc('legend', fontsize=22)
mpl.rcParams['ps.useafm'] = True
mpl.rcParams['pdf.use14corefonts'] = True
mpl.rcParams['text.usetex'] = True
```

# Comparing approaches to subjective probability and statistics

We illustrate various approaches to statistics based on subjective probability. Our test case is the estimation of the mean of a one-dimensional Gaussian. We will demonstrate each approach to subjective probability using one point estimate and one "uncertainty" interval. We will also comment on the mathematical derivation as well as the interpretation of these two outputs for each approach.

The setting is $x_1,\dots,x_N\sim\mathcal{N}(\theta,\sigma^2)$, with $\sigma^2$ known. The statistician's goal is to infer $\theta$, and state his uncertainty about his estimate. 


```python
def generateData(N, mu, sigma2, tag="well-specified"):
    if tag=="well-specified":
        return mu+np.sqrt(sigma2)*npr.randn(N)
```


```python
#npr.seed(1)
thetaTrue = 1 # "true" value of the mean
sigma2True = 1
alpha = .05
gamma = 0.5 # importance of the length of the interval vs. the interval containing the true value
N = 10
x = generateData(N, thetaTrue, sigma2True)
```


```python
plt.hist(x)
plt.show()
```

The dataset is denoted by $\mathcal{D}= \left\lbrace x_1, ..., x_N \right\rbrace$. The corresonding likelihood function is 
$$ \mathcal{L} \left( \theta; \mathcal{D} \right)  = \prod_{i=1}^N \frac{1}{\sigma \sqrt{2\pi}}\exp \left( - \frac{\left( \theta - x_i \right)^2}{2 \sigma^2} \right).$$

The maximum likelihood estimate (MLE) of $\theta$ is the sample mean $\bar{x}$.


```python
intervals = {}
```

## Conjugate Bayes

### The Bayesian paradigm
The Bayesian approach is based on the expected utility principle. 
* First, define a cost function $L(\theta,a)$ that measures how much you suffer from picking action $a$ when the state of the world --here, the variable to infer-- is $\mu$. Actions should be functions from the set of states of the world to a space of "outcomes". For our application, an action should consist in picking an estimate $\hat{\theta}$ and a small enough interval $I$, and the outcomes that matter to know if we've suceeded are $a(\theta) = (\theta-\hat{\theta}, \vert I\vert, 1_{\theta\in I})$. So we need to define $L(\theta,a) = L(\theta,\hat{\theta},I)$. 
* Second, simply choose a pdf $p(\theta)$ over states of the world, your prior, and choose the action $a$ that minimizes the posterior expected loss
$$
\mathbb{E}_\pi L(\theta, a) = \int L(\theta, a) \pi(\theta)d\theta,
$$
where we have denoted by $\pi(\theta)$ the probability distribution function (pdf) proportional to $p(x\vert\theta)p(\theta)$. By Bayes' theorem, this is the conditional pdf of $\theta$ given data $x$, also called the posterior distribution. 

### Choosing a loss function
For the estimate $\hat{\theta}$, we follow here the long tradition of squared losses and we decide to suffer $(\theta-\hat{\theta})^2$ if we output $\hat{\theta}$ while the actual mean of the Gaussian is $\theta$. For the interval, we want to penalize intervals that do not contain the actual $\theta$, while also penalizing intervals that are too wide, to avoid vacuous statements like $I=\mathbb{R}$. This leads us to

$$ L \left (\theta, \hat{\theta}, I \right) = (\theta-\hat{\theta})^2 + \gamma \vert I\vert + 1_{\theta\notin I}. $$

The parameter $\gamma>0$ needs to be chosen by the user, and is interpreted as how much the user favors small intervals over intervals that contain the actual value. There is no need for an additional tradeoff parameter in front of the first term of the right-hand side, since 
$$
\min_{\hat{\theta}, I} \mathbb{E}_\pi L(\theta, \hat{\theta}, I) = \min_{\hat{\theta}}\mathbb{E}_\pi (\theta-\hat{\theta})^2 + \min_I \left[\gamma \vert I \vert + \int 1_{\theta\notin I} \pi(\theta) d\theta \right]. 
$$
Actually this global optimization problem breaks into two simpler ones (wrt to $\theta$ and wrt $I$).

By a property of the mean, the Bayesian estimate $\hat{\theta}$ for our choice of loss is thus the mean of the posterior distribution $\pi$. The Bayesian estimate of $I$ is not obvious at this stage.

### Picking a prior
Now, if our prior $p(\theta)$ is a Gaussian $\mathcal{N}(\mu_\text{prior}, \sigma^2_\text{prior})$, and our likelihood is also Gaussian $\mathcal{N}(\theta,\sigma^2)$, you can check that the posterior $\pi$ is 
\begin{equation}
\mathcal{N}\left( \frac{\frac{\mu_\text{prior}}{\sigma^2_\text{prior}}+\frac{\sum_{i=1}^N x_i}{\sigma^2}}{\frac{1}{\sigma^2_\text{prior}}+\frac{N}{\sigma^2}} ,  \frac{1}{\sigma^2_\text{prior}}+\frac{N}{\sigma^2} \right) .
\label{e:posterior}
\end{equation}
This setting where the posterior has the same functional form (here Gaussian) as the prior is called "conjugate", hence the title of this section. Note that we chose a conjugate prior as a mathematical convenience and because it allows easy interpretation, but being Bayesian does not require it. Note how the mean of the posterior is a convex combination of the prior mean and the average of the data points. `ConjugateBayes.estimate()` simply returns this convex combination. 

Because the posterior is Gaussian, an interval of given width will always contain more posterior mass if it is centered around the posterior mean. Hence, we can limit our search for the best $I$ to intervals $I_\alpha$ centered around the mean, where $\alpha$ indicates that the posterior puts mass $1-\alpha$ on this interval. Computing the endpoints of $I_\alpha$ simply involves calling the `erf` function, and that is what `ConjugateBayes.interval(1-alpha)` does.


```python
import ConjugateBayes as cb
cb = cb.ConjugateBayes(data=x, sigma2Lhd=sigma2True**2, muPrior=10, sigma2Prior=100)
```


```python
# Each approach should be implemented in one class, in a separate python file. Each approach should have the following two methods:
print(cb.estimate()) # yields a point estimate, here the posterior mean
print(cb.interval(1-alpha)) # yields an uncertainty quantification, here a credible interval of posterior mass 1-alpha
```

Now recall that among centered intervals, I want to pick the $I_\alpha$ with the lowest integrated cost 
$$\gamma \vert I_\alpha\vert +  \int \left(1_{\theta\notin I_\alpha}\right) \pi(\theta)d\theta = \gamma \vert I_\alpha\vert + \alpha$$


```python
alphas = np.linspace(1e-5,1-1e-5, 100)
costs_cb = []
for alpha_loc  in alphas:
    # Compute cost for various alphas
    lo, hi = cb.interval(1-alpha_loc)
    costs_cb.append( gamma*(hi-lo) + alpha_loc )

ind = np.argmin(costs_cb) # Find the alpha that minimizes the cost
alphaStar = alphas[ind]
intervals["conjugate Bayes"] = cb.interval(1-alphaStar)
print("The cost is minimized for alpha=", alphaStar)
print("The interval is then", intervals["conjugate Bayes"])
```


```python
# We can visualize the cost as a function of alpha
plt.plot(alphas, costs_cb, color='b') 

# adjust the plotting window and show minimizer
m = np.min(costs_cb) 
M = np.max(costs_cb)
delta = M-m
yMin = m-.5*delta
yMax = M+.5*delta
plt.ylim([yMin,yMax])
plt.vlines(alphaStar,yMin,yMax,linestyle='--',color='blue')

plt.xlabel(r"$\alpha$")
plt.ylabel("cost")
plt.grid()
plt.show()
```

The Bayes action is to report the credible interval with $\alpha$ set to the minizer of this curve! We can also visualize this interval.


```python
for k, res in enumerate(intervals.items()):
    method, interval = res
    plt.plot([k, k], [interval[0], interval[1]], label=method, linewidth=6)
    plt.grid()
    plt.ylabel(r"$\theta$")
    plt.legend()
```

## Finitely additive Bayes

The expected utility principle, along with the rule that says you should update your beliefs by conditioning on observed events, can be axiomatized. One of the most significant work in this direction is known as Savage's axioms, which is usually what statisticians refer to when they claim  that "being Bayesian is being a coherent rational agent". However, Savage's axioms actually only warrant using finitely additive probability measures, and these measures have to be defined for all subsets of the set of states of the world. This is incompatible with using a Gaussian prior over $\theta\in\mathbb{R}$, as there will be nonmeasurable sets for this Gaussian. TBC.

## Likelihood based Belief Function

### Definitions

In the belief function framework, as well as in other ill-known probabilities frameworks, we work with a pair of functions interpreted as lower and upper bounds of the odds of $\theta$ being $\mu$. The lower bound $bel$ is called **belief function** while the upper bound $pl$ is called **plausibility function**. These functions are set functions, which means that their domain is the set of subsets of $\Theta = \mathbb{R}$. For any subset $A\subset \Theta$, the two functions are related by

$$bel \left( A \right) = 1 - pl \left(A^c\right).$$

There are several ways to derive such function. In this notebook, we present belief functions *à la Dempster* which are induced by a random set. Let $(\Omega,\sigma_\Omega,\mu)$ denote a  probability space. Define a mapping $\Gamma$ from $\Omega$ to the set of subsets of $\Theta$. Under measurability conditions, $\Gamma$ is a **random set**. We can define the two following pseudo-inverse for $\Gamma$:
$$ \Gamma^{-1}_{\top}\left(B\right) = \left\lbrace \omega | \Gamma\left(\omega\right)\cap B \neq \emptyset \right\rbrace,\forall B\subseteq \Theta ,\\
\Gamma^{-1}_{\bot}\left(B\right) = \left\lbrace \omega | \Gamma\left(\omega\right)\subseteq B \right\rbrace,\forall B\subseteq \Theta .$$

Then, we have $pl = \mu \circ \Gamma^{-1}_{\top}$ and $bel = \mu \circ \Gamma^{-1}_{\bot} $.

**Example**:

Going back to our inference problem for $\theta$, we can use the likelihood based belief function. The corresponding plausibility function is defined as 
$$ pl \left( \theta \right) = \frac{\mathcal{L} \left( \theta ; \mathcal{D} \right)}{\mathcal{L} \left( \bar{x} ; \mathcal{D} \right)}, \forall \theta \in \mathbb{R} \\
pl \left( H \right) = \underset{\theta \in H}{\sup} pl \left( \theta \right), \forall H \subseteq \mathbb{R}
$$

This model stems from the probability space $(U,\sigma_U,\lambda)$ where $U$ is the unit interval, $\sigma_U$ the Borel sigma-field on this interval and $\lambda$ is the uniform probability measure. Now, by mapping each $u \in U$ to the $\alpha$-cut of function $\frac{\mathcal{L} \left( \theta ; \mathcal{D} \right)}{\mathcal{L} \left( \bar{x} ; \mathcal{D} \right)}$, we span the above plausibility function.

**Interpretation**:
For any event $H$, we can understand the belief and plausibility values as:
* $bel \left( H\right)$ : probability to know that $H$ is true,
* $pl \left( H\right)$ : probability to know that $H^c$ is false,
* $pl \left( H\right) - bel \left( H\right)$ : probability not to know whether $H$ is true.

So we see that the probability mass $pl \left( H\right) - bel \left( H\right)$ is not assigned to any outcome $H$ or $H^c$ in light of the available data.

In the specific case of the Likelihood based belief function, we also have random interval interpretation. In this case, we build the belief function by choosing uniformly a confidence level $\alpha$ and propagate the correponding probability measure to the $\alpha$-cuts of the likelihood function. $bel\left(H\right)$ is then the probability that the random interval is included in $H$ while $pl\left(H\right)$ is the probability that the random interval intersects $H$.

### Belief functions and decision theory

Like the Bayesian paradigm, making decisions with belief functions can also be formalized using expected utility. 
Yet belief functions are non-additive measures and we cannot integrate a cost against them in the sense of Lebesgue integral. For this purpose, we need to resort a non-additive integral now as Choquet integral. Let alone these technicalities, making decisions using a belief function shares the same philosophy as in the Bayesian case:

* First, also pick a loss function. In the sequel, we will go for a different loss function as before but our aim is still the same. We need to pick an estimate of $\theta$ and a small condidence interval which is believed to contain $\theta$. We will also penalize incorrect estimations, large intervals and intervals that miss $\theta$.
* Second, update your prior information with your data. The information drawn from the data is embodied by the likelihood based belief function. The prior information is also a belief function. The two functions are combined through an operation called Dempster's rule which we do not detail here. In the case your prior information is a probability distribution, then Dempster's rule yields a posterior distribution which conincides with the Bayesian approach. In the sequel, we choose a genuinely uninformative prior, i.e. a vacuous belief function ($bel = 1_\Theta$).
* Finally, choose the action that minimizes the integrated cost. However, we have two non additive measures against which it can be integrated :
  $$ (C)\int J\: dbel  \text{ (Conservative approach)},\\
  (C)\int J\: dpl  \text{ (Optimistic approach)},$$
  
  where $(C)\int$ is the Choquet integral symbol.
  
Just like Bayesians do with Savage's representation theorem, BF practionners invoke Gilboa's representation theorem to claim that the BF approach is compliant with rational decision making. As proved by Gilboa, appraising the cost of an action based on expected utility relying on non-additive measures translates into a preference relation among actions with axioms similar to Savage's. Two important restrictions appears in both Savage's and Gilboa's axiomatizations: the utility function must be bounded (we will comply to this requirement in the sequel) and the measure must have convex range.

### Choosing a loss function

The non-additivity of belief function leads to non-intuitive calculations for people who are familiar with standard probability calculus. In particular, we can use the $0-1$ loss $L_{0-1}$ for the estimation of $\theta$, while this would be pointless within Lebesgue integrals. In the likelihood belief function case, we have 

# John, can you justify here separating the cost minimization into the point estimate and the interval estimate? I think you had an argument based on comonotonicity

$$ (C)\int L_{0-1}\: dbel = 1 - pl \left( \hat{\theta} \right) ,\\
  (C)\int L_{0-1}\: dpl = 1.$$

In the belief case, we obtain $\hat{\theta} = \bar{x}$, while plausibility cannot draw any conclusion. We thus use the belief function for the estimation of $\theta$, and the method  `BaselineBF.estimate` returns the MLE estimate.

For the interval estimation problem, we will use the same cost function as before, and we obtain
$$ (C)\int \vert I\vert + 1_{\theta\notin I} \: dbel = \gamma \vert I \vert ,\\
  (C)\int  \vert I\vert + 1_{\theta\notin I}\: dpl = \gamma \vert I \vert + \sup_{\theta\notin I}pl(\theta).$$
  
The integral wrt the belief function is always minimized when $\vert I\vert=0$, leaving $I$ undetermined. For the plausibility function, the likelihood being Gaussian, it is enough to minimize the expected cost among intervals that are upper level sets of the likelihood, that is 
$$ I_\alpha = \{\theta/ pl(\theta)<=\alpha\}.$$
We minimize over $\alpha$ below.

To conclude, note that we cannot use the same non-additive measure for both estimation problems. This is in contrast with the Bayesian approach.


```python
import BaselineBF as bf
bf = bf.BaselineBF(data=x, sigma2Lhd=sigma2True)
```


```python
print(bf.estimate()) # yields a point estimate, here the MLE
```

  Now, we are still interested in the same interval cost function $J \left(\alpha \right)$ 
  
  The method `BaselineBF.interval` returns the estimated "confidence" interval.


```python
print(bf.interval(1-alpha))
```


```python
costs_cons = []
costs_optm = []
for alpha_loc  in alphas:
    lo, hi = bf.interval(1-alpha_loc)
    costs_cons.append( gamma*(hi-lo))
    costs_optm.append( gamma*(hi-lo) + alpha_loc )

yMin = 0
yMax = np.max(costs_optm)
plt.ylim([yMin,yMax])

# Conservative approach
ind = np.argmin(costs_cons) # Find the alpha that minimizes the cost
alphaStar = alphas[ind]
plt.vlines(alphaStar,yMin,yMax,linestyle='--',color='red')
print("The conservative cost is minimized for alpha=", alphaStar)
lo, hi = bf.interval(1-alphaStar)
intervals["conservative likelihood-based BF"] = [lo, hi]
print("The interval is then", intervals["conservative likelihood-based BF"])

# Optimistic approach
ind = np.argmin(costs_optm) # Find the alpha that minimizes the cost
alphaStar = alphas[ind]
plt.vlines(alphaStar,yMin,yMax,linestyle='--',color='green')
print("The optimistic cost is minimized for alpha=", alphaStar)
lo, hi = bf.interval(1-alphaStar)
intervals["optimistic likelihood-based BF"] = [lo, hi]
print("The interval is then", intervals["optimistic likelihood-based BF"])

# plot the costs
plt.plot(alphas, costs_cons, label="conservative", color='r')
plt.plot(alphas, costs_optm, label="optimistic", color='g')
plt.plot(alphas, costs_cb, label="conjugate Bayes", color='b')

# show the minimzers

plt.xlabel(r"$\alpha$")
plt.ylabel("cost")
plt.legend()
plt.show()
```

We focus on the optimistic cost as the conservative does not provide a sufficient level of information. In this experiment, the optimistic cost is minimized for the same value of $\alpha$ as in the conjugate Bayes case. This may be purely conincidental as other values of $\gamma$ may not lead to such a result. Obviously, when $\gamma$ is small, both costs become closer. 

We also see that the optimistic cost curve is not convex while the conjugate Bayes curve seems to be.

## Comparing intervals


```python
for k, res in enumerate(intervals.items()):
    method, interval = res
    plt.plot([k, k], [interval[0], interval[1]], label=method, linewidth=6)
    plt.grid()
    plt.legend()
```


```python
print(intervals)
```


```python

```


```python

```

060087bd7079883f1402ad53f9f231f99fef4ff1
ipynb
Jupyter Notebook
tutorial_notebook/.ipynb_checkpoints/tutorial-checkpoint.ipynb
rbardenet/subjective-fits
e0115ba9f527b62158846fe68b1e8a9b86573d42
MIT
2018-03-11T11:35:32.000Z
2018-04-27T08:26:39.000Z
tutorial_notebook/.ipynb_checkpoints/tutorial-checkpoint.ipynb
rbardenet/subjective-fits
e0115ba9f527b62158846fe68b1e8a9b86573d42
MIT
tutorial_notebook/.ipynb_checkpoints/tutorial-checkpoint.ipynb
rbardenet/subjective-fits
e0115ba9f527b62158846fe68b1e8a9b86573d42
MIT
Qwen/Qwen-72B
1. YES
2. YES