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()