# lampproef.py - is de nieuwe lamp echt beter? Analyse op het niveau van de proefeenheid (het vakje)
import csv
import math
from collections import defaultdict
from scipy import stats

planten = defaultdict(list)                          # (vakje, lamp) -> versgewichten van de zes planten
with open("lampproef.csv", newline="") as f:
    for rij in csv.DictReader(f):
        planten[(int(rij["vakje"]), rij["lamp"])].append(float(rij["versgewicht_g"]))


def gem(x):
    return sum(x) / len(x)


def var(x):                                          # steekproefvariantie (delen door n - 1)
    m = gem(x)
    return sum((v - m) ** 2 for v in x) / (len(x) - 1)


def welch(a, b):
    va, vb = var(a) / len(a), var(b) / len(b)
    se = math.sqrt(va + vb)
    t = (gem(a) - gem(b)) / se
    nu = (va + vb) ** 2 / (va ** 2 / (len(a) - 1) + vb ** 2 / (len(b) - 1))
    return gem(a) - gem(b), se, t, nu, 2 * stats.t.sf(abs(t), nu)


vak = {lamp: [gem(g) for (v, l), g in sorted(planten.items()) if l == lamp] for lamp in ("oud", "nieuw")}
for lamp in ("oud", "nieuw"):
    print("%-5s vakjes: %s  gem %.1f g, s %.1f g" % (lamp, " ".join("%.1f" % x for x in vak[lamp]),
                                                 gem(vak[lamp]), math.sqrt(var(vak[lamp]))))

d, se, t, nu, p = welch(vak["nieuw"], vak["oud"])
tk = stats.t.ppf(0.975, nu)
print("vakjes (n = 6 + 6): verschil %.1f g, se %.2f g, t = %.2f, nu = %.1f, p = %.3f" % (d, se, t, nu, p))
print("   95 %%-interval: %.1f tot %.1f g   (scipy: p = %.3f)"
      % (d - tk * se, d + tk * se, stats.ttest_ind(vak["nieuw"], vak["oud"], equal_var=False).pvalue))

# zo NIET: elke plant als onafhankelijke herhaling tellen (pseudoreplicatie)
alle = {lamp: [x for (v, l), g in planten.items() if l == lamp for x in g] for lamp in ("oud", "nieuw")}
d2, se2, t2, nu2, p2 = welch(alle["nieuw"], alle["oud"])
print("planten (n = 36 + 36): verschil %.1f g, se %.2f g, t = %.2f, p = %.4f  <- te zelfzeker" % (d2, se2, t2, p2))