# Part 1: meta-analytic r (Macnamara et al. 2014) as r^2. Part 2: synthetic population
# (illustrative weights) in which practice helps EVERY person, yet explains little BETWEEN them.
import math, random
from statistics import correlation, fmean, pvariance

print("Part 1  meta-analytic correlations (Macnamara et al. 2014)")
for dom, r in (("overall", .35), ("games", .51), ("music", .46), ("sports", .42),
               ("education", .21), ("professions", .05)):
    print(f"   {dom:11s} r = {r:.2f} -> r^2 = {r * r:5.1%} explained, {1 - r * r:5.1%} not")

rng = random.Random(23)
people = []
for _ in range(5000):
    hours = math.exp(rng.gauss(math.log(3000), 0.6))       # accumulated practice, log-normal
    quality = rng.gauss(0, 1)                                # how well the hours are used
    start = rng.gauss(0, 1)                                  # earlier start (standardised)
    ability = rng.gauss(0, 1)                                # prior ability, aptitude
    other = rng.gauss(0, 1)                                  # teaching, health, luck, ...
    gain = 1.0 + 0.4 * quality                               # each person's return per log-hour
    perf = gain * (math.log(hours) - math.log(300)) + 0.35 * start + 0.8 * ability + 0.9 * other
    people.append((math.log(hours), quality, start, ability, other, perf, gain))

perf = [p[5] for p in people]
def r2(i):
    return correlation([p[i] for p in people], perf) ** 2

print("\nPart 2  synthetic population of 5,000 learners (illustrative)")
print(f"   within a person: +1 log-unit of practice (x2.7 hours) adds on average "
      f"{fmean(p[6] for p in people):.2f} performance units")
sd = math.sqrt(pvariance(perf))
print(f"   between people: SD of performance {sd:.2f} units")
print(f"   so going from 300 to 3,000 hours moves a typical person {math.log(10) / sd:.2f} SD")
for name, i in (("log practice hours", 0), ("practice quality", 1), ("starting age", 2),
                ("prior ability", 3), ("everything else", 4)):
    print(f"   share explained by {name:18s} alone: {r2(i):5.1%}")
qp = [p[0] * p[6] for p in people]                           # hours weighted by quality
print(f"   share explained by quality-weighted hours:  "
      f"{correlation(qp, perf) ** 2:5.1%}")
# Part 1  meta-analytic correlations (Macnamara et al. 2014)
#    overall     r = 0.35 -> r^2 = 12.2% explained, 87.8% not
#    games       r = 0.51 -> r^2 = 26.0% explained, 74.0% not
#    music       r = 0.46 -> r^2 = 21.2% explained, 78.8% not
#    sports      r = 0.42 -> r^2 = 17.6% explained, 82.4% not
#    education   r = 0.21 -> r^2 =  4.4% explained, 95.6% not
#    professions r = 0.05 -> r^2 =  0.3% explained, 99.8% not
#
# Part 2  synthetic population of 5,000 learners (illustrative)
#    within a person: +1 log-unit of practice (x2.7 hours) adds on average 1.00 performance units
#    between people: SD of performance 1.69 units
#    so going from 300 to 3,000 hours moves a typical person 1.36 SD
#    share explained by log practice hours alone: 12.1%
#    share explained by practice quality   alone: 29.6%
#    share explained by starting age       alone:  4.1%
#    share explained by prior ability      alone: 22.6%
#    share explained by everything else    alone: 28.5%
#    share explained by quality-weighted hours:  37.0%

# One true relation (log practice vs performance, r = 0.45, synthetic data) seen through
# (1) range restriction: only the top performers are studied; (2) measurement: practice
# estimated with large error, or recalled in a way that tracks present performance.
import math, random
from statistics import correlation

rng = random.Random(7)
pop = [(lp, 0.45 * lp + math.sqrt(1 - 0.45 ** 2) * rng.gauss(0, 1))
       for lp in (rng.gauss(0, 1) for _ in range(20000))]
r2 = lambda rows: correlation([a for a, _ in rows], [b for _, b in rows]) ** 2
print(f"whole population: practice explains {r2(pop):.1%}")
ranked = sorted(pop, key=lambda r: r[1])
for top in (0.10, 0.02):
    print(f"   top {top:.0%} performers only: {r2(ranked[int(len(ranked) * (1 - top)):]):.1%}")
for label, err, bias in (("large recall error", 1.0, 0.0), ("recall tracks performance", 1.0, 0.5)):
    obs = [(lp + err * rng.gauss(0, 1) + bias * pf, pf) for lp, pf in pop]
    print(f"   {label}: {r2(obs):.1%}")
# whole population: practice explains 20.0%
#    top 10% performers only: 3.4%
#    top 2% performers only: 1.3%
#    large recall error: 9.8%
#    recall tracks performance: 33.8%