# Which generations are DIRECTLY exposed, and which is the first that is not?
# The exposed individual is F0; germ cells are exposed with the body carrying
# them. Part 2 is a toy model with ILLUSTRATIVE numbers: a germline mark must
# escape two waves of reprogramming per generation, whereas a behavioural
# channel (maternal care -> pup -> her own care) is rebuilt every generation.
EXPOSED = {"adult male": ["F0", "F1 (germ cells)"],
           "non-pregnant female": ["F0", "F1 (germ cells)"],
           "pregnant female": ["F0", "F1 (fetus)", "F2 (germ cells)"]}
for case, gens in EXPOSED.items():
    print(f"{case:19s}: exposed {', '.join(gens)}; first unexposed F{len(gens)}")

esc = 0.10 * 0.50          # escape both waves (illustrative)
fidelity = 0.70            # behavioural channel: share kept per generation
print("share of an induced signal left (illustrative)")
print(" generation | germline mark, escape 0.10 x 0.50 | behavioural channel 0.70")
for g in range(1, 5):
    print(f"     F{g}     | {esc ** g:33.6f} | {fidelity ** g:24.3f}")
# adult male         : exposed F0, F1 (germ cells); first unexposed F2
# non-pregnant female: exposed F0, F1 (germ cells); first unexposed F2
# pregnant female    : exposed F0, F1 (fetus), F2 (germ cells); first unexposed F3
# share of an induced signal left (illustrative)
#  generation | germline mark, escape 0.10 x 0.50 | behavioural channel 0.70
#      F1     |                          0.050000 |                    0.700
#      F2     |                          0.002500 |                    0.490
#      F3     |                          0.000125 |                    0.343
#      F4     |                          0.000006 |                    0.240

# A simulated cross-fostering experiment in the style of Francis et al. (1999).
# Daughter's adult licking-and-grooming score (LG, standardised) =
#   0.6 x LG of her REARING mother + noise   (behavioural channel), or
#   0.6 x LG of her BIRTH mother + noise     (germline channel).
# Half the pups are swapped to an unrelated mother. Numbers illustrative.
import random, statistics as st
random.seed(22)

def experiment(channel, n=400, k=0.6):
    bio = [random.gauss(0, 1) for _ in range(n)]
    rear = list(bio)
    idx = list(range(n // 2, n))                # second half: swapped mothers
    for i, j in zip(idx, idx[1:] + idx[:1]): rear[i] = bio[j]
    src = rear if channel == "behavioural" else bio
    dau = [k * s + random.gauss(0, (1 - k * k) ** 0.5) for s in src]
    f = slice(n // 2, n)
    return st.correlation(bio[f], dau[f]), st.correlation(rear[f], dau[f])

print("cross-fostered daughters only (200 per hypothesis)")
print(" hypothesis  | r with biological mother | r with foster mother")
for ch in ("behavioural", "germline"):
    rb, rr = experiment(ch)
    print(f" {ch:11s} | {rb:24.2f} | {rr:20.2f}")
# cross-fostered daughters only (200 per hypothesis)
#  hypothesis  | r with biological mother | r with foster mother
#  behavioural |                     0.07 |                 0.63
#  germline    |                     0.65 |                -0.07

# Why small human methylation studies need replication. Simulated studies
# compare "exposed offspring" (n = 20) with controls (n = 10) at 20 CpG sites,
# using Welch's t-test. Case 1: no true difference anywhere. Case 2: one site
# has a true difference of 0.3 SD. Sizes and effect are illustrative.
import math, random, statistics as st
random.seed(7)

def betacf(a, b, x):            # continued fraction, incomplete beta function
    qab, qap, qam, c, d = a + b, a + 1, a - 1, 1.0, 1 - (a + b) * x / (a + 1)
    d = 1 / d; h = d
    for m in range(1, 200):
        m2 = 2 * m
        for aa in (m * (b - m) * x / ((qam + m2) * (a + m2)),
                   -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))):
            d = 1 / (1 + aa * d); c = 1 + aa / c; h *= d * c
    return h

def p_two_sided(t, df):         # P(|T| > |t|), Student's t with df d.o.f.
    x = df / (df + t * t)
    lb = math.lgamma(df / 2 + 0.5) - math.lgamma(df / 2) - math.lgamma(0.5)
    bt = math.exp(lb + (df / 2) * math.log(x) + 0.5 * math.log(1 - x))
    if x < (df / 2 + 1) / (df / 2 + 1.5):
        return bt * betacf(df / 2, 0.5, x) / (df / 2)
    return 1 - bt * betacf(0.5, df / 2, 1 - x) / 0.5

def welch(x, y):
    vx, vy = st.variance(x) / len(x), st.variance(y) / len(y)
    t = (st.fmean(x) - st.fmean(y)) / math.sqrt(vx + vy)
    df = (vx + vy) ** 2 / (vx ** 2 / (len(x) - 1) + vy ** 2 / (len(y) - 1))
    sp = math.sqrt((st.variance(x) + st.variance(y)) / 2)
    return p_two_sided(t, df), (st.fmean(x) - st.fmean(y)) / sp

def study(true_d, n1=20, n0=10, sites=20):
    res = []
    for s in range(sites):
        shift = true_d if s == 0 else 0.0
        res.append(welch([random.gauss(shift, 1) for _ in range(n1)],
                         [random.gauss(0, 1) for _ in range(n0)]))
    return res

R = 1000
null = [study(0.0) for _ in range(R)]
share = sum(any(p < .05 for p, _ in s) for s in null) / R
print(f"no true effect: studies with at least one site p < .05: {share:.0%}")
true = [study(0.3) for _ in range(R)]
hits = [d for s in true for (p, d) in s[:1] if p < .05]
print(f"true 0.3 SD at one site: power at that site {len(hits) / R:.0%};")
print(f"  mean estimated effect when significant {st.fmean(hits):.2f} SD")
big = sum(study(0.3, 200, 100, 1)[0][0] < .05 for _ in range(500)) / 500
print(f"with n = 200 vs 100 the same site is detected in {big:.0%} of studies")
# no true effect: studies with at least one site p < .05: 65%
# true 0.3 SD at one site: power at that site 12%;
#   mean estimated effect when significant 0.91 SD
# with n = 200 vs 100 the same site is detected in 66% of studies