import numpy as np, matplotlib.pyplot as plt, warnings
warnings.filterwarnings("ignore")
import pymc as pm, arviz as az
from scipy import stats
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"

cross=np.array([23.5,12,21,22,19.125,21.5,22.125,20.375,18.25,21.625,23.25,21,22.125,23,12])
selff=np.array([17.375,20.375,20,20,18.375,18.625,18.625,15.25,16.5,18,16.25,18,12.75,15.5,18])
d = cross - selff; n = len(d)

t_stat, t_p = stats.ttest_rel(cross, selff)
se = d.std(ddof=1)/np.sqrt(n); tcrit = stats.t.ppf(0.975, n-1)
ci = (d.mean()-tcrit*se, d.mean()+tcrit*se)
print(f"{n} pairs; mean difference = {d.mean():.3f} in, sd = {d.std(ddof=1):.3f}, SE = {se:.3f}")
print(f"paired t-test: t = {t_stat:.4f}, two-sided p = {t_p:.4f}")
print(f"95% t-interval = [{ci[0]:.3f}, {ci[1]:.3f}]  -- clears zero by {ci[0]:.3f} inches")
print(f"\nthe two negative pairs are {d[d<0]} against {(d>0).sum()} positives -- the outliers the t-test must absorb")

# the frequentist benchmark this page argues with, recomputed here rather than quoted
import itertools
signs = np.array(list(itertools.product([1,-1], repeat=n)))
perm  = (signs*d).mean(1)
n_ge  = int((perm >= d.mean()).sum())
p_one, p_two = n_ge/len(perm), float((np.abs(perm) >= abs(d.mean())).mean())
print(f"\nFisher's exact permutation test, all {2**n:,} sign patterns enumerated:")
print(f"  one-sided p = {p_one:.5f}   ({n_ge:,}/{2**n:,} patterns give a mean this large)")
print(f"  two-sided p = {p_two:.5f}")
print(f"  paired t-test, two-sided p = {t_p:.5f}")
print()
print("Compared like with like, the two disagree on the conventional threshold: the assumption-free")
print(f"permutation test gives {p_two:.4f} and FAILS at 0.05, while the normal-theory t-test gives {t_p:.4f}")
print("and passes. Three thousandths separate them. Neither is wrong -- the result is sitting on the")
print("boundary, which is the whole reason this page exists.")

g++ not available, if using conda: `conda install gxx`

15 pairs; mean difference = 2.617 in, sd = 4.718, SE = 1.218
paired t-test: t = 2.1480, two-sided p = 0.0497
95% t-interval = [0.004, 5.229]  -- clears zero by 0.004 inches

the two negative pairs are [-8.375 -6.   ] against 13 positives -- the outliers the t-test must absorb

Fisher's exact permutation test, all 32,768 sign patterns enumerated:
  one-sided p = 0.02634   (863/32,768 patterns give a mean this large)
  two-sided p = 0.05267
  paired t-test, two-sided p = 0.04970

Compared like with like, the two disagree on the conventional threshold: the assumption-free
permutation test gives 0.0527 and FAILS at 0.05, while the normal-theory t-test gives 0.0497
and passes. Three thousandths separate them. Neither is wrong -- the result is sitting on the
boundary, which is the whole reason this page exists.

with pm.Model() as m_normal:
    tau   = pm.Normal("tau", 0, 10)
    sigma = pm.HalfNormal("sigma", 10)
    pm.Normal("obs", tau, sigma, observed=d)
    idata = pm.sample(2000, tune=1000, chains=4, cores=1, progressbar=False, random_seed=0)

post = idata.posterior["tau"].values.ravel()
lo, hi = np.percentile(post, [2.5, 97.5])
print(f"posterior mean tau = {post.mean():.3f} in   (t-test point estimate {d.mean():.3f})")
print(f"95% credible interval = [{lo:.3f}, {hi:.3f}]   (t-interval [{ci[0]:.3f}, {ci[1]:.3f}])")
print(f"P(tau > 0 | data)     = {(post>0).mean():.4f}")
print(f"P(tau > 1 inch)       = {(post>1).mean():.4f}")
print(f"max R-hat = {float(az.summary(idata)['r_hat'].max()):.4f}")
print()
print("Note what happened to the verdict. The t-interval EXCLUDES zero by 0.004 inches; this credible")
print(f"interval INCLUDES it, running from {lo:.3f}. And P(tau>0)={(post>0).mean():.4f} sits just below the 0.975 that a")
print("one-sided 0.025 test would need. The same data, a prior nobody would call informative -- sd 10 on an")
print("effect of 2.6 inches -- and the significance evaporates. That is not a Bayesian trick; it is what a")
print("result sitting four thousandths of an inch from the boundary was always going to do under any")
print("reanalysis. The posterior also answers the question actually asked: P(tau>0) is the probability the")
print("effect is positive, where 0.0497 is the probability of data this extreme if it were exactly zero.")

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (4 chains in 1 job)

NUTS: [tau, sigma]

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 1 seconds.

posterior mean tau = 2.573 in   (t-test point estimate 2.617)
95% credible interval = [-0.050, 5.256]   (t-interval [0.004, 5.229])
P(tau > 0 | data)     = 0.9721
P(tau > 1 inch)       = 0.8826
max R-hat = 1.0000

Note what happened to the verdict. The t-interval EXCLUDES zero by 0.004 inches; this credible
interval INCLUDES it, running from -0.050. And P(tau>0)=0.9721 sits just below the 0.975 that a
one-sided 0.025 test would need. The same data, a prior nobody would call informative -- sd 10 on an
effect of 2.6 inches -- and the significance evaporates. That is not a Bayesian trick; it is what a
result sitting four thousandths of an inch from the boundary was always going to do under any
reanalysis. The posterior also answers the question actually asked: P(tau>0) is the probability the
effect is positive, where 0.0497 is the probability of data this extreme if it were exactly zero.

scales = [0.5, 1, 2, 3, 5, 10, 25, 100]
rows=[]
for s in scales:
    with pm.Model():
        tau   = pm.Normal("tau", 0, s)
        sigma = pm.HalfNormal("sigma", 10)
        pm.Normal("obs", tau, sigma, observed=d)
        it = pm.sample(2000, tune=1000, chains=2, cores=1, progressbar=False, random_seed=1)
    p = it.posterior["tau"].values.ravel()
    rows.append((s, p.mean(), np.percentile(p,2.5), np.percentile(p,97.5), (p>0).mean()))

print(f"  {'prior sd':>9} {'post mean':>10} {'2.5%':>8} {'97.5%':>8} {'P(tau>0)':>10}")
for s,mn,l,h,pp in rows:
    print(f"  {s:>9} {mn:>10.3f} {l:>8.3f} {h:>8.3f} {pp:>10.4f}")

pp = np.array([r[4] for r in rows]); mns = np.array([r[1] for r in rows])
fig,ax=plt.subplots(1,2,figsize=(13,4.2))
ax[0].semilogx(scales, pp, "o-", color=BLUE, lw=2)
ax[0].axhline(0.975, color=GREEN, ls="--", lw=1, label="0.975  (one-sided 0.025)")
ax[0].axhline(1-t_p/2, color=RED, ls=":", lw=1.5, label=f"1 - p/2 = {1-t_p/2:.3f}  (t-test)")
ax[0].set_xlabel("prior sd on tau (inches)"); ax[0].set_ylabel("P(tau > 0 | data)")
ax[0].set_title("How much of the verdict is the prior?"); ax[0].legend(fontsize=8); ax[0].set_ylim(0.5,1.0)
ax[1].semilogx(scales, mns, "o-", color=PURP, lw=2, label="posterior mean")
ax[1].fill_between(scales, [r[2] for r in rows], [r[3] for r in rows], color=PURP, alpha=.15, label="95% credible")
ax[1].axhline(d.mean(), color=GREY, ls="--", lw=1, label=f"sample mean {d.mean():.2f}")
ax[1].axhline(0, color="k", lw=.6)
ax[1].set_xlabel("prior sd on tau (inches)"); ax[1].set_ylabel("effect (inches)")
ax[1].set_title("Shrinkage toward the sceptical prior"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()

print(f"At the flat end (s=100) P(tau>0) = {pp[-1]:.4f}, essentially the one-sided t-test answer of {1-t_p/2:.4f}.")
print(f"At s=1 -- a prior saying the effect is probably under an inch -- it is {pp[1]:.4f}, and the point")
print(f"estimate has shrunk from {d.mean():.2f} to {mns[1]:.2f} inches.")
print("Read the curve rather than a single number. P(tau>0) holds above 0.94 for any prior sd of 2 or more --")
print("that is, any prior that does not rule out an effect the size of the one observed. Below that it falls")
print(f"away quickly: {pp[1]:.3f} at sd 1 and {pp[0]:.3f} at sd 0.5, where the prior asserts the effect is almost")
print("certainly under half an inch. So the direction is robust to ignorance but not to active scepticism,")
print("and the magnitude is fragile throughout -- the point estimate runs from 0.31 to 2.62 across the sweep.")
print("Fifteen pots cannot settle how big the effect is. They can only say it is probably positive.")

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (2 chains in 1 job)

NUTS: [tau, sigma]

Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 1 seconds.

We recommend running at least 4 chains for robust computation of convergence diagnostics

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (2 chains in 1 job)

NUTS: [tau, sigma]

Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 1 seconds.

We recommend running at least 4 chains for robust computation of convergence diagnostics

with pm.Model() as m_robust:
    tau   = pm.Normal("tau", 0, 10)
    sigma = pm.HalfNormal("sigma", 10)
    nu    = pm.Gamma("nu", alpha=2, beta=0.1)
    pm.StudentT("obs", nu=nu, mu=tau, sigma=sigma, observed=d)
    id_r = pm.sample(2000, tune=1000, chains=4, cores=1, progressbar=False, random_seed=2)

pr = id_r.posterior["tau"].values.ravel(); nur = id_r.posterior["nu"].values.ravel()
lo_r, hi_r = np.percentile(pr,[2.5,97.5])
print(f"  {'model':22s} {'effect':>8} {'95% interval':>22} {'P(tau>0)':>10}")
print(f"  {'paired t-test':22s} {d.mean():>8.3f} {'[%.3f, %.3f]'%ci:>22} {1-t_p/2:>10.4f}")
print(f"  {'Bayes, normal':22s} {post.mean():>8.3f} {'[%.3f, %.3f]'%(lo,hi):>22} {(post>0).mean():>10.4f}")
print(f"  {'Bayes, Student-t':22s} {pr.mean():>8.3f} {'[%.3f, %.3f]'%(lo_r,hi_r):>22} {(pr>0).mean():>10.4f}")
print(f"\nposterior median nu = {np.median(nur):.1f}  (low = heavy tails; above ~30 is effectively normal)")
print(f"P(nu < 10 | data) = {(nur<10).mean():.3f}")

fig,ax=plt.subplots(1,2,figsize=(13,4.2))
ax[0].hist(post,bins=60,density=True,color=BLUE,alpha=.55,label="normal likelihood")
ax[0].hist(pr,bins=60,density=True,color=ORANGE,alpha=.55,label="Student-t likelihood")
ax[0].axvline(0,color="k",lw=.8); ax[0].axvline(d.mean(),color=GREY,ls="--",lw=1,label="sample mean")
ax[0].set_xlabel("effect of cross-fertilization (inches)"); ax[0].set_ylabel("posterior density")
ax[0].set_title("The outliers, handled two ways"); ax[0].legend(fontsize=8)
ax[1].hist(nur,bins=60,color=PURP,alpha=.7); ax[1].axvline(np.median(nur),color=RED,lw=1.5,label=f"median {np.median(nur):.1f}")
ax[1].set_xlabel("nu (degrees of freedom)"); ax[1].set_ylabel("posterior draws")
ax[1].set_title("How heavy do the tails need to be?"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()

print("\nThe robust model moves the estimate UP, not down, and tightens the interval: down-weighting the two")
print("negative pairs removes variance the t-test had to carry. Fisher avoided the normality assumption by")
print("refusing to make one; the Bayesian route makes it and then lets the data say how wrong it is.")

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (4 chains in 1 job)

NUTS: [tau, sigma, nu]

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 2 seconds.

  model                    effect           95% interval   P(tau>0)
  paired t-test             2.617         [0.004, 5.229]     0.9751
  Bayes, normal             2.573        [-0.050, 5.256]     0.9721
  Bayes, Student-t          2.855         [0.183, 5.374]     0.9806

posterior median nu = 15.2  (low = heavy tails; above ~30 is effectively normal)
P(nu < 10 | data) = 0.308

The robust model moves the estimate UP, not down, and tightens the interval: down-weighting the two
negative pairs removes variance the t-test had to carry. Fisher avoided the normality assumption by
refusing to make one; the Bayesian route makes it and then lets the data say how wrong it is.

Causal Inference II(a) — Bayesian Analysis of a Small Experiment¶

Darwin's maize, where 15 pairs leave the prior doing visible work¶

1. The data and the frequentist benchmark¶

2. The posterior, and what it answers that a $p$-value does not¶

3. The prior-sensitivity curve — the thing the frequentist analysis cannot draw¶

4. Robustness — the assumption Fisher was avoiding¶

5. Summary¶