# What small hormone-behaviour effects mean in practice.
import math
from statistics import NormalDist
N = NormalDist()

def from_d(d):
    """Standardised mean difference -> chance that a random member of group 1 scores higher than a
    random member of group 2 (common-language effect size), and the overlap of the two distributions."""
    cl = N.cdf(d / 2 ** 0.5)
    overlap = 2 * N.cdf(-abs(d) / 2)
    return cl, overlap

print("d      P(random member of higher group is higher)   overlap of distributions")
for d in (0.1, 0.2, 0.5, 0.8):
    cl, ov = from_d(d)
    print(f"{d:4.1f}   {100 * cl:5.1f} %                                     {100 * ov:5.1f} %")

r = -0.061          # testosterone x cortisol interaction, Dekkers et al. 2019 meta-analysis
print(f"\nr = {r}: explains {100 * r * r:.2f} % of the variance in status-relevant behaviour")
# People needed to detect such an effect once, two-sided alpha 0.05, 80 % power (Fisher z approximation)
z = 0.5 * math.log((1 + abs(r)) / (1 - abs(r)))
n = ((N.inv_cdf(0.975) + N.inv_cdf(0.80)) / z) ** 2 + 3
print(f"participants needed for 80 % power to detect r = {abs(r)}: {math.ceil(n)}")
# d      P(random member of higher group is higher)   overlap of distributions
#  0.1    52.8 %                                      96.0 %
#  0.2    55.6 %                                      92.0 %
#  0.5    63.8 %                                      80.3 %
#  0.8    71.4 %                                      68.9 %
#
# r = -0.061: explains 0.37 % of the variance in status-relevant behaviour
# participants needed for 80 % power to detect r = 0.061: 2108

# Toy model of two coupled axes during a 60-minute contest (arbitrary units, rest = 1).
# Cortisol C rises with appraised demand; testosterone T production is boosted by a contest the person
# appraises as a winnable challenge and is inhibited by high cortisol (HPA -> HPG cross-talk, Viau 2002).
# Parameters are illustrative; effect sizes in real people are much smaller and more variable (see text).
def contest(c_drive, t_boost, minutes=240, dt=0.5):
    C = T = 1.0
    hC, hT, KI = 1 / 60, 1 / 45, 2.0          # clearance per minute; inhibition constant
    out, t = {}, 0.0
    while t <= minutes:
        on = t < 60
        dC = hC * ((1 + (c_drive if on else 0)) - C)
        # production: boosted during a challenge, inhibited by cortisol (equals 1 at C = 1 without a boost)
        prod = (1 + (t_boost if on else 0)) * (1 + (1 / KI) ** 2) / (1 + (C / KI) ** 2)
        dT = hT * (prod - T)
        C += dt * dC; T += dt * dT
        for m in (30, 60, 120, 240):
            if abs(t - m) < dt / 2:
                out[m] = (C, T)
        t += dt
    return out

print(f"{'appraisal':32s}  min     C     T   T/C")
for name, cd, tb in (("challenge (demands <= resources)", 0.8, 0.4),
                     ("threat (demands > resources)", 2.5, 0.0),
                     ("no contest", 0.0, 0.0)):
    o = contest(cd, tb)
    for m in (60, 120, 240):
        C, T = o[m]
        print(f"{name if m == 60 else '':32s} {m:4d}: {C:5.2f} {T:5.2f} {T / C:5.2f}")
# appraisal                         min     C     T   T/C
# challenge (demands <= resources)   60:  1.50  1.15  0.76
#                                   120:  1.18  0.95  0.81
#                                   240:  1.02  0.98  0.95
# threat (demands > resources)       60:  2.57  0.71  0.28
#                                   120:  1.58  0.68  0.43
#                                   240:  1.08  0.92  0.85
# no contest                         60:  1.00  1.00  1.00
#                                   120:  1.00  1.00  1.00
#                                   240:  1.00  1.00  1.00