|
2 | 2 | Bayesian Estimation Supersedes the T-Test
|
3 | 3 |
|
4 | 4 | This model replicates the example used in:
|
5 |
| -Kruschke, John. (2012) Bayesian estimation supersedes the t test. Journal of Experimental Psychology: General. |
| 5 | +Kruschke, John. (2012) Bayesian estimation supersedes the t test. |
| 6 | +Journal of Experimental Psychology: General. |
6 | 7 |
|
7 |
| -The original pymc2 implementation was written by Andrew Straw and can be found here: https://github.com/strawlab/best |
| 8 | +The original pymc2 implementation was written by Andrew Straw |
| 9 | +and can be found here: |
| 10 | +https://github.com/strawlab/best |
8 | 11 |
|
9 | 12 | Ported to PyMC3 by Thomas Wiecki (c) 2015.
|
10 | 13 | """
|
11 | 14 |
|
12 | 15 | import numpy as np
|
| 16 | + |
13 | 17 | import pymc3 as pm
|
14 | 18 |
|
15 |
| -drug = (101,100,102,104,102,97,105,105,98,101,100,123,105,103,100,95,102,106, |
16 |
| - 109,102,82,102,100,102,102,101,102,102,103,103,97,97,103,101,97,104, |
17 |
| - 96,103,124,101,101,100,101,101,104,100,101) |
18 |
| -placebo = (99,101,100,101,102,100,97,101,104,101,102,102,100,105,88,101,100, |
19 |
| - 104,100,100,100,101,102,103,97,101,101,100,101,99,101,100,100, |
20 |
| - 101,100,99,101,100,102,99,100,99) |
| 19 | +drug = (101, 100, 102, 104, 102, 97, 105, 105, 98, 101, 100, |
| 20 | + 123, 105, 103, 100, 95, 102, 106, |
| 21 | + 109, 102, 82, 102, 100, 102, 102, 101, 102, 102, 103, |
| 22 | + 103, 97, 97, 103, 101, 97, 104, |
| 23 | + 96, 103, 124, 101, 101, 100, 101, 101, 104, 100, 101) |
| 24 | +placebo = (99, 101, 100, 101, 102, 100, 97, 101, 104, 101, 102, |
| 25 | + 102, 100, 105, 88, 101, 100, 104, 100, 100, 100, |
| 26 | + 101, 102, 103, 97, 101, 101, |
| 27 | + 100, 101, 99, 101, 100, 100, |
| 28 | + 101, 100, 99, 101, 100, 102, 99, 100, 99) |
21 | 29 |
|
22 | 30 | y1 = np.array(drug)
|
23 | 31 | y2 = np.array(placebo)
|
24 | 32 | y = np.concatenate((y1, y2))
|
25 | 33 |
|
26 |
| -mu_m = np.mean( y ) |
27 |
| -mu_p = 0.000001 * 1/np.std(y)**2 |
| 34 | +mu_m = np.mean(y) |
| 35 | +mu_p = 0.000001 * 1 / np.std(y) ** 2 |
28 | 36 |
|
29 |
| -sigma_low = np.std(y)/1000 |
30 |
| -sigma_high = np.std(y)*1000 |
| 37 | +sigma_low = np.std(y) / 1000 |
| 38 | +sigma_high = np.std(y) * 1000 |
31 | 39 |
|
32 | 40 | with pm.Model() as model:
|
33 |
| - group1_mean = pm.Normal('group1_mean', mu=mu_m, tau=mu_p, testval=y1.mean()) |
34 |
| - group2_mean = pm.Normal('group2_mean', mu=mu_m, tau=mu_p, testval=y2.mean()) |
35 |
| - group1_std = pm.Uniform('group1_std', lower=sigma_low, upper=sigma_high, testval=y1.std()) |
36 |
| - group2_std = pm.Uniform('group2_std', lower=sigma_low, upper=sigma_high, testval=y2.std()) |
37 |
| - nu = pm.Exponential('nu_minus_one', 1/29.) + 1 |
38 |
| - |
39 |
| - lam1 = group1_std**-2 |
40 |
| - lam2 = group2_std**-2 |
41 |
| - |
42 |
| - group1 = pm.StudentT('drug', nu=nu, mu=group1_mean, lam=lam1, observed=y1) |
43 |
| - group2 = pm.StudentT('placebo', nu=nu, mu=group2_mean, lam=lam2, observed=y2) |
44 |
| - |
45 |
| - diff_of_means = pm.Deterministic('difference of means', group1_mean - group2_mean) |
46 |
| - diff_of_stds = pm.Deterministic('difference of stds', group1_std - group2_std) |
47 |
| - effect_size = pm.Deterministic('effect size', diff_of_means / pm.sqrt((group1_std**2 + group2_std**2) / 2)) |
| 41 | + group1_mean = pm.Normal('group1_mean', mu=mu_m, tau=mu_p, |
| 42 | + testval=y1.mean()) |
| 43 | + group2_mean = pm.Normal('group2_mean', mu=mu_m, tau=mu_p, |
| 44 | + testval=y2.mean()) |
| 45 | + group1_std = pm.Uniform('group1_std', lower=sigma_low, upper=sigma_high, |
| 46 | + testval=y1.std()) |
| 47 | + group2_std = pm.Uniform('group2_std', lower=sigma_low, upper=sigma_high, |
| 48 | + testval=y2.std()) |
| 49 | + nu = pm.Exponential('nu_minus_one', 1 / 29.) + 1 |
| 50 | + |
| 51 | + lam1 = group1_std ** -2 |
| 52 | + lam2 = group2_std ** -2 |
| 53 | + |
| 54 | + group1 = pm.StudentT( |
| 55 | + 'drug', nu=nu, mu=group1_mean, lam=lam1, observed=y1) |
| 56 | + group2 = pm.StudentT( |
| 57 | + 'placebo', nu=nu, mu=group2_mean, lam=lam2, |
| 58 | + observed=y2) |
| 59 | + |
| 60 | + diff_of_means = pm.Deterministic( |
| 61 | + 'difference of means', group1_mean - |
| 62 | + group2_mean) |
| 63 | + diff_of_stds = pm.Deterministic( |
| 64 | + 'difference of stds', |
| 65 | + group1_std - group2_std) |
| 66 | + effect_size = pm.Deterministic( |
| 67 | + 'effect size', |
| 68 | + diff_of_means / pm.sqrt( |
| 69 | + (group1_std ** |
| 70 | + 2 + group2_std**2) / 2)) |
48 | 71 |
|
49 | 72 | step = pm.NUTS()
|
50 | 73 |
|
| 74 | + |
51 | 75 | def run(n=3000):
|
52 | 76 | if n == "short":
|
53 | 77 | n = 500
|
54 | 78 | with model:
|
55 | 79 | trace = pm.sample(n, step)
|
56 | 80 |
|
57 |
| - burn = n/10 |
| 81 | + burn = n / 10 |
58 | 82 |
|
59 |
| - pm.traceplot(trace[burn:]); |
| 83 | + pm.traceplot(trace[burn:]) |
60 | 84 | pm.plots.summary(trace[burn:])
|
61 | 85 |
|
| 86 | + |
62 | 87 | if __name__ == '__main__':
|
63 | 88 | run()
|
0 commit comments