# Bayesian causal inference for a flash and a beep (Koerding et al. 2007). C = 1: one
# source; C = 2: two. Group-fitted values from the paper: sd_V 2.14, sd_A 9.2 deg, prior
# on position N(0, 12.3^2), p_common 0.28. Beep at 0 (prior centre), flash at +D. Bias of
# the heard position = mean estimate / D, in % of the way to the flash.
import math, random
SV, SA, SP, PC = 2.14, 9.2, 12.3, 0.28

def post_common(xv, xa, sv=SV, sa=SA, sp=SP, pc=PC):
    v1, va, vp = sv * sv, sa * sa, sp * sp
    k = v1 * va + v1 * vp + va * vp
    l1 = math.exp(-((xv - xa) ** 2 * vp + xv ** 2 * va + xa ** 2 * v1) / (2 * k)) \
        / (2 * math.pi * math.sqrt(k))
    l2 = math.exp(-xv ** 2 / (2 * (v1 + vp)) - xa ** 2 / (2 * (va + vp))) \
        / (2 * math.pi * math.sqrt((v1 + vp) * (va + vp)))
    return pc * l1 / (pc * l1 + (1 - pc) * l2)

def heard(xv, xa, pc=PC):
    p1 = post_common(xv, xa, pc=pc)
    fused = (xv / SV ** 2 + xa / SA ** 2) / (1 / SV ** 2 + 1 / SA ** 2 + 1 / SP ** 2)
    alone = (xa / SA ** 2) / (1 / SA ** 2 + 1 / SP ** 2)
    return p1 * fused + (1 - p1) * alone, p1

if __name__ == "__main__":
    rng = random.Random(2007)
    print("D (deg)  P(one source)   bias of heard beep:  causal   forced   none"
          "   causal, p_c 0.8")
    for d in (0, 5, 10, 15, 20, 30, 40):
        res = [(rng.gauss(d, SV), rng.gauss(0.0, SA)) for _ in range(20000)]
        p1 = sum(post_common(xv, xa) for xv, xa in res) / len(res)
        if d == 0:
            print(f"{d:6.0f} {p1:12.2f}          (no conflict)")
            continue
        b = [round(100 * sum(heard(xv, xa, pc)[0] for xv, xa in res) / len(res) / d)
             for pc in (PC, 1.0, 0.0, 0.8)]
        print(f"{d:6.0f} {p1:12.2f} {b[0]:23d}% {b[1]:7d}% {b[2]:5d}% {b[3]:11d}%")
# D (deg)  P(one source)   bias of heard beep:  causal   forced   none   causal, p_c 0.8
#      0         0.32          (no conflict)
#      5         0.31                      20%      93%     1%          62%
#     10         0.26                      17%      92%     0%          55%
#     15         0.20                      12%      92%     0%          44%
#     20         0.14                       9%      92%     0%          33%
#     30         0.05                       3%      92%     0%          13%
#     40         0.01                       1%      92%     0%           3%

# How much of "the" McGurk effect is the person and how much the clip? Illustrative
# measurement model: P(fusion) = logistic(person + day + clip). Person ~ N(0, 2.5) and
# day-to-day wobble ~ N(0, 0.7) on the logit scale (both ASSUMED); 12 clips tuned so the
# population fusion rate runs from 17 % to 58 %, the range reported by Basu Mallick et
# al. (2015). 10 trials per clip per session.
import math, random
from statistics import correlation
logistic = lambda z: 1 / (1 + math.exp(-z))

def tune(target, people):
    lo, hi = -10.0, 10.0
    for _ in range(60):
        mid = (lo + hi) / 2
        rate = sum(logistic(p + mid) for p in people) / len(people)
        lo, hi = (mid, hi) if rate < target else (lo, mid)
    return hi

if __name__ == "__main__":
    rng = random.Random(1976)
    people = [rng.gauss(0, 2.5) for _ in range(2000)]
    clips = [tune(0.17 + (0.58 - 0.17) * j / 11, people) for j in range(12)]
    def session():
        out = []
        for p in people:
            q = p + rng.gauss(0, 0.7)
            out.append([sum(rng.random() < logistic(q + c) for _ in range(10)) / 10
                        for c in clips])
        return out
    s1, s2 = session(), session()
    each = [sum(r) / 12 for r in s1]
    rates = [sum(r[j] for r in s1) / len(s1) for j in range(12)]
    print("clip fusion rates (population):", " ".join(f"{x:.0%}" for x in rates))
    print(f"people, all 12 clips: {sum(e < 0.1 for e in each) / len(each):.0%} fuse on"
          f" <10 % of trials, {sum(e > 0.9 for e in each) / len(each):.0%} on >90 %,"
          f" range {min(each):.0%} to {max(each):.0%}")
    for name, j in (("weakest clip", 0), ("strongest clip", 11)):
        share = sum(r[j] >= 0.5 for r in s1) / len(s1)
        print(f"'has the McGurk effect' (>= 5 of 10 fusions) on the {name}: {share:.0%}")
    r = correlation(each, [sum(x) / 12 for x in s2])
    r1 = correlation([x[0] for x in s1], [x[0] for x in s2])
    print(f"test-retest r: mean over 12 clips {r:.2f}; one weak clip alone {r1:.2f}")
# clip fusion rates (population): 18% 21% 25% 29% 32% 36% 40% 43% 47% 51% 55% 58%
# people, all 12 clips: 27% fuse on <10 % of trials, 7% on >90 %, range 0% to 100%
# 'has the McGurk effect' (>= 5 of 10 fusions) on the weakest clip: 16%
# 'has the McGurk effect' (>= 5 of 10 fusions) on the strongest clip: 62%
# test-retest r: mean over 12 clips 0.91; one weak clip alone 0.78