import numpy as np, matplotlib.pyplot as plt, pymc as pm, seaborn as sns
%matplotlib inline

def model():
    sample_0 = [8.254828927, 8.485694524, 10.58058423, 5.221356325]
    sample_1 = [7.921107724, 9.744301571, 5.750874082, 8.421883012, 11.09813068]
    pooled_data = np.hstack((sample_0, sample_1))

    # hyper-parameters
    S = 1000 * np.std(pooled_data)
    M = np.mean(pooled_data)

    L = .001 * np.std(pooled_data)
    H = 1000 * np.std(pooled_data)

    # prior
    mu = pm.Normal('mu_1', M, S**-2, value=[M,M])
    sigma = pm.Uniform('sigma_1', L, H, value=[np.std(pooled_data), np.std(pooled_data)])

    # likelihood
    obs_1 = pm.Normal('obs_1', mu=mu[0], tau=sigma[0]**-2, value=sample_0, observed=True)
    obs_2 = pm.Normal('obs_2', mu=mu[1], tau=sigma[1]**-2, value=sample_1, observed=True)

    # derived quantities of interest
    d_mu = mu[0] - mu[1]
    d_sigma = sigma[0] - sigma[1]
    @pm.deterministic
    def effect_size(d_mu=d_mu, sigma=sigma):
        return d_mu / np.sqrt(.5*(sigma[0]**2 + sigma[1]**2))
    @pm.deterministic
    def outside_rope(effect_size=effect_size):
        return -.1 < effect_size < .1
    
    return locals()

vars = model()

m = pm.MCMC(vars)
m.use_step_method(pm.AdaptiveMetropolis, [m.mu, m.sigma])
m.sample(100000, 50000, 10)

pm.Matplot.plot(m)

# estimate of probability that effect size
# is outside region of practical equivalence
m.outside_rope.trace().mean()  # not "significant"