Implementation of 8 schools, with a twist.

hierarchical_normal.py

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import arviz as az


N = 50_000
treatment_conversions = np.array([541, 557, 559, 556, 530, 532, 516, 532, 528, 544, 519, 552])
control_conversions = np.array([496, 524, 486, 500, 516, 475, 507, 475, 490, 506, 512, 489])
log_estimated_rr = np.log(treatment_conversions / control_conversions)
sd_log_estimated_rr = 1/treatment_conversions + 1/control_conversions - 2/N
experiment = ["Experiment_{i}" for i in range(treatment_conversions.size)]

with pm.Model(coords={'experiment':experiment}) as model:
    
    mu = pm.Normal('mu',mu=0, sigma=1)
    tau = pm.HalfCauchy('tau', 2.5)
    z = pm.Normal('z',0, 1, dims='experiment')
    
    theta = pm.Deterministic('theta', mu + tau * z)
    
    log_RR = pm.Normal('log_RR', theta, sd_log_estimated_rr, observed=log_estimated_rr)
    
with model:
    trace = pm.sample(cores=4, chains=4)