# Detection bias in official records (the question of Widom et al. 2015).
# 100,000 parents, 15 % with a documented maltreatment history. TRUE rates of
# maltreatment of their children 12 % vs 9 % (ILLUSTRATIVE). Records catch a
# true case with probability 0.60 if the family is already known to child
# protection, 0.30 otherwise; a confidential survey catches 70 % in both
# groups, with 2 % false positives.
import random
rng = random.Random(1967)

def odds_ratio(a, n1, b, n0):
    return (a / (n1 - a)) / (b / (n0 - b))

N, share = 100000, 0.15
hist = [rng.random() < share for _ in range(N)]
true = [rng.random() < (0.12 if h else 0.09) for h in hist]
rec = [t and rng.random() < (0.60 if h else 0.30) for h, t in zip(hist, true)]
rep = [(t and rng.random() < 0.70) or (not t and rng.random() < 0.02)
       for t in true]
n1 = sum(hist); n0 = N - n1
for name, v in (("true maltreatment", true), ("official records", rec),
                ("self-report", rep)):
    a = sum(x for x, h in zip(v, hist) if h)
    b = sum(x for x, h in zip(v, hist) if not h)
    print(f"{name:18s}: {a / n1:5.1%} vs {b / n0:5.1%}   "
          f"odds ratio {odds_ratio(a, n1, b, n0):.2f}")
print("records alone would suggest a 'cycle' twice as strong as the true one")
# true maltreatment : 11.7% vs  9.0%   odds ratio 1.34
# official records  :  7.0% vs  2.7%   odds ratio 2.74
# self-report       :  9.9% vs  8.0%   odds ratio 1.26
# records alone would suggest a 'cycle' twice as strong as the true one

# Two-state Markov chain for harsh parenting (H) across generations.
# P(H | parent H) = 0.30 (Kaufman and Zigler's rough "about 30 %", as an
# example); P(H | parent N) = 0.10 (ILLUSTRATIVE). A buffer (partner,
# mentor, programme) reaches a share phi of adults raised in H families and
# removes a share k = 0.7 (ILLUSTRATIVE) of the excess risk.

def chain(phi, k=0.7, p_hh=0.30, p_nh=0.10, gens=4):
    eff = (1 - phi) * p_hh + phi * (p_nh + (p_hh - p_nh) * (1 - k))
    lam = eff - p_nh; x, path = 1.0, []
    for _ in range(gens):
        x = x * eff + (1 - x) * p_nh; path.append(x)
    return eff, lam, path, p_nh / (1 - lam)

print("lineage that starts with a harsh parent (generation 0)")
print(" buffer | P(H|H) |  lam  | P(H) in generations 1-4    | long-run P(H)")
for phi in (0.0, 0.3, 0.6, 0.9):
    eff, lam, path, ss = chain(phi)
    print(f"   {phi:.1f}  |  {eff:.3f} | {lam:.3f} |"
          f" {', '.join(f'{p:.3f}' for p in path)} |    {ss:.3f}")
print("with no continuity at all (P(H|H) = P(H|N) = 0.10): long-run P(H) = "
      f"{chain(0.0, p_hh=0.10)[3]:.3f}")
print("even with no buffer, 70 % of children of harsh parents do not become")
print("harsh parents")
# lineage that starts with a harsh parent (generation 0)
#  buffer | P(H|H) |  lam  | P(H) in generations 1-4    | long-run P(H)
#    0.0  |  0.300 | 0.200 | 0.300, 0.160, 0.132, 0.126 |    0.125
#    0.3  |  0.258 | 0.158 | 0.258, 0.141, 0.122, 0.119 |    0.119
#    0.6  |  0.216 | 0.116 | 0.216, 0.125, 0.115, 0.113 |    0.113
#    0.9  |  0.174 | 0.074 | 0.174, 0.113, 0.108, 0.108 |    0.108
# with no continuity at all (P(H|H) = P(H|N) = 0.10): long-run P(H) = 0.100
# even with no buffer, 70 % of children of harsh parents do not become
# harsh parents