# Reparametizing Gelman's rat tumour example
**URL:** https://discourse.pymc.io/t/reparametizing-gelmans-rat-tumour-example/1417
**Category:** Questions
**Created:** [June 23, 2018, 9:18pm UTC](https://discourse.pymc.io/t/reparametizing-gelmans-rat-tumour-example/1417 "2018-06-23T21:18:58Z")
**Posts on this page:** 8
**Page:** 1
### Author:  [@dpananos](https://discourse.pymc.io/u/dpananos)
#### Post date: [June 23, 2018, 9:18pm UTC](https://discourse.pymc.io/t/reparametizing-gelmans-rat-tumour-example/1417/1 "2018-06-23T21:18:58Z")
This question concerns Gelman’s Rat Tumor example in BDA 3 (bottom of page 109).
I’m trying to replicate the inference made in this example. Here is some code to get started:
```
import numpy as np
import pandas as pd
import pymc3 as pm
import theano.tensor as T
import matplotlib.pyplot as plt
# rat data (BDA3, p. 102)
n_tumour = np.array([
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1,
1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 5, 2,
5, 3, 2, 7, 7, 3, 3, 2, 9, 10, 4, 4, 4, 4, 4, 4, 4,
10, 4, 4, 4, 5, 11, 12, 5, 5, 6, 5, 6, 6, 6, 6, 16, 15,
15, 9, 4
])
n_rat = np.array([
20, 20, 20, 20, 20, 20, 20, 19, 19, 19, 19, 18, 18, 17, 20, 20, 20,
20, 19, 19, 18, 18, 25, 24, 23, 20, 20, 20, 20, 20, 20, 10, 49, 19,
46, 27, 17, 49, 47, 20, 20, 13, 48, 50, 20, 20, 20, 20, 20, 20, 20,
48, 19, 19, 19, 22, 46, 49, 20, 20, 23, 19, 22, 20, 20, 20, 52, 46,
47, 24, 14
])
with pm.Model() as model:
#Uninformative prior for alpha and beta
u = pm.Uniform('u', 0, 1)
v = pm.Uniform('v', 0,1)
alpha = pm.Deterministic('a',u/T.pow(v,2))
beta = pm.Deterministic('b',(1-u)/T.pow(v,2))
theta = pm.Beta('theta', alpha=alpha, beta = beta)
p = pm.Binomial('y',p=theta,observed = n_tumour, n=n_rat)
trace = pm.sample(50_000, cores = 2, step = pm.NUTS())
```
Gelman choses a uniform hyperprior on (\dfrac{a}{a+b}, (a+b)^{-1/2}), so let u = \dfrac{a}{a+b} and v = (a+b)^{-1/2} and solve appropriately.
I keep getting the error
```
There were 4377 divergences after tuning. Increase `target_accept` or reparameterize.
There were 2637 divergences after tuning. Increase `target_accept` or reparameterize.
The number of effective samples is smaller than 10% for some parameters.
```
so I was wondering how best to reparameterize the problem so that the chains converge. Am I reparameterizing incorrectly by using `pm.Deterministic`?
---
### Author:  [@junpenglao](https://discourse.pymc.io/u/junpenglao)
#### Post date: [June 24, 2018, 8:39am UTC](https://discourse.pymc.io/t/reparametizing-gelmans-rat-tumour-example/1417/2 "2018-06-24T08:39:38Z")
Hmm, if I understand correctly what BDA described is uniform on the domain of the parameter, not just on [0, 1]
---
### Author:  [@junpenglao](https://discourse.pymc.io/u/junpenglao)
#### Post date: [June 24, 2018, 2:39pm UTC](https://discourse.pymc.io/t/reparametizing-gelmans-rat-tumour-example/1417/3 "2018-06-24T14:39:11Z")
Actually, the problem is because you didnt model a vector of `theta`:
```python
with pm.Model() as model:
# Uninformative prior for alpha and beta
u = pm.Uniform('u', 0, 1)
v = pm.Uniform('v', 0, 1)
alpha = pm.Deterministic('a', u / tt.pow(v, 2))
beta = pm.Deterministic('b', (1 - u) / tt.pow(v, 2))
# theta_i with i=1,2,...,J
theta = pm.Beta('theta', alpha=alpha, beta=beta, shape=len(n_rat))
p = pm.Binomial('y', p=theta, observed=n_tumour, n=n_rat)
trace = pm.sample(1000, tune=1000)
pm.traceplot(trace);
```
Alternatively, following the equation (5.9) in BDA 3 you can directly specify the prior of \alpha, \beta using a `DensityDist`:
```python
import pymc3.distributions.transforms as tr
def logp_ab(value):
return tt.log(tt.pow(tt.sum(value), -5/2))
with pm.Model() as model:
# Uninformative prior for alpha and beta
ab = pm.DensityDist('ab', logp_ab, transform=tr.log,
shape=2,
testval=np.asarray([1., 1.]))
theta = pm.Beta('theta', alpha=ab[0], beta=ab[1], shape=len(n_rat))
p = pm.Binomial('y', p=theta, observed=n_tumour, n=n_rat)
```
---
### Author:  [@dpananos](https://discourse.pymc.io/u/dpananos)
#### Post date: [June 24, 2018, 3:44pm UTC](https://discourse.pymc.io/t/reparametizing-gelmans-rat-tumour-example/1417/4 "2018-06-24T15:44:52Z")
Thanks!
---
### Author:  [@dpananos](https://discourse.pymc.io/u/dpananos)
#### Post date: [June 25, 2018, 10:06pm UTC](https://discourse.pymc.io/t/reparametizing-gelmans-rat-tumour-example/1417/5 "2018-06-25T22:06:06Z")
When I run the model with `transform = tr.log`, I underestimate the value of alpha and beta.
According to Gelman, the expected means for alpha and beta should be 2.4 and 14.3 respectively. I can verify this by directly computing the posterior.
When I run this model, I estimate the means for alpha and beta to be 1.8 and 10.8. Running without `transform = tr.log` yields estimate closer to the true values. Since you take the log in `logp_ab`, is `transform = tr.log` needed again?
---
### Author:  [@junpenglao](https://discourse.pymc.io/u/junpenglao)
#### Post date: [June 26, 2018, 5:26am UTC](https://discourse.pymc.io/t/reparametizing-gelmans-rat-tumour-example/1417/6 "2018-06-26T05:26:22Z")
If you don’t apply the transformation, NUTS will give you divergences when the purpose value went out of bound of the support. I guess that’s why in BDA3 they have another equation for unbound parameters log(a/b) and log(a+b) (eq 5.10)
Here the transformation added a jacobian which biased the estimation. To get one with no transformation you can first create the free parameter using `HalfFlat` and then add the logp using potential:
```python
def logp_ab(value):
return tt.log(tt.pow(tt.sum(value), -5/2))
with pm.Model() as model:
# Uninformative prior for alpha and beta
ab = pm.HalfFlat('ab',
shape=2,
testval=np.asarray([1., 1.]))
pm.Potential('p(a, b)', logp_ab(ab))
theta = pm.Beta('theta', alpha=ab[0], beta=ab[1], shape=len(n_rat))
p = pm.Binomial('y', p=theta, observed=n_tumour, n=n_rat)
trace = pm.sample(1000, tune=2000, nuts_kwargs={'target_accept': .95})
```
---
### Author:  [@junpenglao](https://discourse.pymc.io/u/junpenglao)
#### Post date: [June 26, 2018, 5:27am UTC](https://discourse.pymc.io/t/reparametizing-gelmans-rat-tumour-example/1417/7 "2018-06-26T05:27:16Z")
For your PR that adds this example, could you also add a paragraph about this?
---
### Author:  [@dpananos](https://discourse.pymc.io/u/dpananos)
#### Post date: [June 26, 2018, 1:37pm UTC](https://discourse.pymc.io/t/reparametizing-gelmans-rat-tumour-example/1417/8 "2018-06-26T13:37:15Z")
Yea, I’ll absolutely add a paragraph about this. Thanks for all your help, I’ll list you as an author on the example.