## This notebook is meant to capture experience with Jeffery's Priors
https://en.wikipedia.org/wiki/Jeffreys_prior




```python
import pymc3 as pm
import theano
import numpy as np
import scipy.stats as ss
import seaborn as sns

sns.set(font_scale=1.5)
```

# Guassian mean

$f(x\mid\mu) = \frac{e^{-(x - \mu)^2 / 2\sigma^2}}{\sqrt{2 \pi \sigma^2}}$

with $\sigma$ fixed, the Jeffreys prior for the mean $\mu$ is
: $\begin{align} p(\mu) & \propto \sqrt{I(\mu)}
= \sqrt{\operatorname{E}\!\left[ \left( \frac{d}{d\mu} \log f(x\mid\mu) \right)^2\right]}
= \sqrt{\operatorname{E}\!\left[ \left( \frac{x - \mu}{\sigma^2} \right)^2 \right]} \\\\
& = \sqrt{\int_{-\infty}^{+\infty} f(x\mid\mu) \left(\frac{x-\mu}{\sigma^2}\right)^2 dx}
=  \sqrt{1/\sigma^2} \propto 1.\end{align}

# Gaussian sd

$f(x\mid\sigma) = \frac{e^{-(x - \mu)^2 / 2 \sigma^2}}{\sqrt{2 \pi \sigma^2}}$,

with $\mu$ fixed, the Jeffreys prior for the standard deviation $\sigma > 0$ is
: $\begin{align}p(\sigma) & \propto \sqrt{I(\sigma)}
= \sqrt{\operatorname{E}\!\left[ \left( \frac{d}{d\sigma} \log f(x\mid\sigma) \right)^2\right]}
= \sqrt{\operatorname{E}\!\left[ \left( \frac{(x - \mu)^2-\sigma^2}{\sigma^3} \right)^2 \right]} \\\\
& = \sqrt{\int_{-\infty}^{+\infty} f(x\mid\sigma)\left(\frac{(x-\mu)^2-\sigma^2}{\sigma^3}\right)^2 dx}
= \sqrt{\frac{2}{\sigma^2}}
\propto \frac{1}{\sigma}.
\end{align} 

### Example model


```python
data = np.random.normal(10, 2, size=20)
sns.distplot(data)
```


```python
def jeffreys(value):
    return -theano.tensor.log(value)

with pm.Model() as model:
    sigma = pm.DensityDist('sigma', jeffreys, testval=1)
    mu = pm.DensityDist('mu', jeffreys, testval=1)
    

```

# Poisson


```python
data = np.random.poisson(20, size=10)
sns.distplot(data)
```


```python
def jeffreys(value):
    return -theano.tensor.log(value)

with pm.Model() as model:
    mean = pm.DensityDist('mean', jeffreys, testval=1)
    poi = pm.Poisson('Poisson', mean, observed=data)
    
    mean_flat = pm.Uniform('mean_flat', 0, 100)
    poi_flat = pm.Poisson('Poisson_flat', mean_flat, observed=data)

    trace = pm.sample(10000, tune=1000)
    
```

    Auto-assigning NUTS sampler...
    Initializing NUTS using jitter+adapt_diag...
    Multiprocess sampling (4 chains in 4 jobs)
    NUTS: [mean_flat, mean]
    Sampling 4 chains, 0 divergences: 100%|██████████| 44000/44000 [00:08<00:00, 5069.74draws/s]



```python
pm.summary(trace)
```




<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>mean</th>
      <th>sd</th>
      <th>hpd_3%</th>
      <th>hpd_97%</th>
      <th>mcse_mean</th>
      <th>mcse_sd</th>
      <th>ess_mean</th>
      <th>ess_sd</th>
      <th>ess_bulk</th>
      <th>ess_tail</th>
      <th>r_hat</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>mean</th>
      <td>18.900</td>
      <td>1.377</td>
      <td>16.334</td>
      <td>21.474</td>
      <td>0.007</td>
      <td>0.005</td>
      <td>35207.0</td>
      <td>34612.0</td>
      <td>35565.0</td>
      <td>27445.0</td>
      <td>1.0</td>
    </tr>
    <tr>
      <th>mean_flat</th>
      <td>19.015</td>
      <td>1.389</td>
      <td>16.434</td>
      <td>21.613</td>
      <td>0.007</td>
      <td>0.005</td>
      <td>34386.0</td>
      <td>34217.0</td>
      <td>34498.0</td>
      <td>28372.0</td>
      <td>1.0</td>
    </tr>
  </tbody>
</table>
</div>




```python
pm.traceplot(trace)
```


```python
pm.plot_posterior(trace)
```


```python

```


```python

```


```python

```


```python

```


```python

```


```python

```


```python

```


```python

```


```python

```


```python
# Make custom Jeffreys prior, returns log prior
def jeffreys(value):
    return -theano.tensor.log(value)

# Instantiate model
model = pm.Model()
with model:
    # Priors for unknown model parameters
    sigma = pm.DensityDist('sigma', jeffreys, testval=1)
    slope = pm.Uniform('slope', lower=-1, upper=1, testval=0.5)
    intercept = pm.Uniform('intercept', lower=-10, upper=10, testval=0)

    # Expected value of outcome
    mu = slope * df.parent_bd + intercept

    # Likelihood (sampling distribution) of observations
    y_obs = pm.Normal('y_obs', mu=mu, sd=sigma, observed=df.offspring_bd)

    # draw 10000 posterior samples for burn-in
    trace = pm.sample(10000)


```


```python
pm.summary(trace, ['slope', 'intercept', 'sigma'])
```


```python
pm.traceplot(trace, ['slope', 'intercept', 'sigma']);
```


```python

```


```python

```


```python

```


```python

```


```python

```


```python

```

8e3fdc0ec51b8415651f7d4d2f3de7ab94cfce4a
ipynb
Jupyter Notebook
Jeffreys prior/Jeffreys prior examples.ipynb
balarsen/pymc_learning
e4a077d492af6604a433433e64b835ce4ed0333a
BSD-3-Clause
Jeffreys prior/Jeffreys prior examples.ipynb
balarsen/pymc_learning
e4a077d492af6604a433433e64b835ce4ed0333a
BSD-3-Clause
Jeffreys prior/Jeffreys prior examples.ipynb
balarsen/pymc_learning
e4a077d492af6604a433433e64b835ce4ed0333a
BSD-3-Clause
2017-05-23T16:38:55.000Z
2017-05-23T16:38:55.000Z
Qwen/Qwen-72B
1. YES
2. YES