<jupyter_start><jupyter_text>First we have to get the bangladesh dataset. [Feather](https://github.com/wesm/feather) is perfect for this. Thanks to <NAME> and <NAME> for producing this invaluable set of tools.The code to export from R should look something like this, assuming you've installed rethinking and feather packages and their dependencies:```Rlibrary(rethinking)library (feather)data(bangladesh)path <- "."write_feather("bangladesh.feather", path)```That's all the R we'll need for now.Next step is to read in the data<jupyter_code># Read the bangledesh data from home directory
bangladesh = ft.read_dataframe("bangladesh.feather")
bangladesh.head()<jupyter_output><empty_output><jupyter_text>From here on we're following Richard McElreath's blog.We set the random seed so we can reproduce the results as required and we sample 100 records from the bangladesh set.<jupyter_code>seed = np.random.RandomState(seed=1)

b = bangladesh.sample(100, random_state=seed)
b.head()<jupyter_output><empty_output><jupyter_text>We only need contraception use tied to district. The data dataframe will hold this information.<jupyter_code>dat = pd.DataFrame()
dat['y'] = b['use.contraception']
dat['district_id'] = pd.Categorical(b['district']).codes 
dat.head()<jupyter_output><empty_output><jupyter_text>Our first model is the **centered parameterization**:$$\begin{align}y_i & \sim \operatorname{Bernoulli}(p_i)\\\operatorname{logit}(p_i) & = \alpha_{unit[i]}\\\alpha_j & \sim \operatorname{Normal}(a, \tau)\\a & \sim \operatorname{Normal}(0, 10)\\\tau & \sim \operatorname{Half-Cauchy}(0, 1)\\\end{align}$$<jupyter_code>with pm.Model() as c_model:
    tau = pm.HalfCauchy('tau',1)
    a = pm.Normal('a', 0, 10)
    alpha = pm.Normal('alpha', a, tau, shape=len(b['district'].unique()))
    p = pm.math.invlogit(alpha[dat['district_id'].values])
    y = pm.Bernoulli('y', p, observed=dat['y'])
    
#     for RV in c_model.basic_RVs:
#         print(RV.name, RV.logp(c_model.test_point))
    trace_c = pm.sample(1000, tune=1000)<jupyter_output>Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
/home/paulc/anaconda3/envs/py36/lib/python3.6/site-packages/pymc3/model.py:384: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if not np.issubdtype(var.dtype, float):
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [alpha, a, tau_log__]
 57%|█████▋    | 1135/2000 [00:05<00:03, 219.45it/s]INFO (theano.gof.compilelock): Waiting for existing lock by process '16081' (I am process '16082')
INFO (theano.gof.compilelock): To manually release the lock, delete /home/paulc/.theano/compiledir_Linux-4.13--generic-x86_64-with-debian-stretch-sid-x86_64-3.6.3-64/lock_dir
 59%|█████▊    | 1172/2000 [00:05<00:03, 222.26it/s]INFO (theano.gof.compilelock): Waiting for existing lock by process '16081' (I am process '16083')
INFO (theano.gof.compilelock): To manually release the lock, delete /home/p[...]<jupyter_text>As with the R Stan version, there's a whole bunch of divergences in this log. This suggests that the results may not be reliable. Our second model is the **non-centered parameterization**:$$\begin{align}y_i & \sim \operatorname{Bernoulli}(p_i)\\\operatorname{logit}(p_i) & = a + \tau z_{unit[i]}\\z_j & \sim \operatorname{Normal}(0, 1)\\a & \sim \operatorname{Normal}(0, 10)\\\tau & \sim \operatorname{Half-Cauchy}(0, 1)\\\end{align}$$<jupyter_code>with pm.Model() as nc_model:
    # pymc3 HalfCauchy only uses the second (scale) parameter
    tau = pm.HalfCauchy('tau', beta=1)
    a = pm.Normal('a', 0, 10)
    z = pm.Normal('z', 0, 1, shape=len(b['district'].unique()))
    p = pm.math.invlogit(a + tau * z[dat['district_id'].values])
    y = pm.Bernoulli('y', p, observed=dat['y'])
    
#     for RV in c_model.basic_RVs:
#         print(RV.name, RV.logp(c_model.test_point))
    trace_nc = pm.sample(1000, tune=1000)<jupyter_output>Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
/home/paulc/anaconda3/envs/py36/lib/python3.6/site-packages/pymc3/model.py:384: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if not np.issubdtype(var.dtype, float):
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [z, a, tau_log__]
100%|██████████| 2000/2000 [00:06<00:00, 302.29it/s]
There were 3 divergences after tuning. Increase `target_accept` or reparameterize.
The acceptance probability does not match the target. It is 0.7201418506887862, but should be close to 0.8. Try to increase the number of tuning steps.
There were 3 divergences after tuning. Increase `target_accept` or reparameterize.
The acceptance probability does not match the target. It is 0.6577402127331399, but should be close to 0.8. Try to increase the number of tuning steps.<jupyter_text>This log looks much better. That suggests already that we might see more reliable results from this approach. Let's take a look at the traceplots:<jupyter_code>_ = pm.plots.traceplot(trace_c)
d = pm.diagnostics.effective_n(trace_c)
print(f"n_eff('a') = {d['a']}\nn_eff('tau') = {d['tau']}")
_ = pm.plots.traceplot(trace_nc)
d = pm.diagnostics.effective_n(trace_nc)
print(f"n_eff('a') = {d['a']}\nn_eff('tau') = {d['tau']}")<jupyter_output>n_eff('a') = 3319.0
n_eff('tau') = 1321.0
{'fraction_non_alphanumeric': 0.14754398826979473, 'fraction_numerical': 0.041788856304985335, 'mean_word_length': 5.0768374164810695, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 24, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '172225', 'license': 'MIT', 'n_tokens_mistral': 2135, 'n_tokens_neox': 1914, 'n_words': 509, 'path': 'Centered vs Non Centered - McElreath.ipynb', 'repo_name': 'paulc00/mlreparam', 'size': '5470'}