# 验证: 共线特征下 "预测稳定" 与 "归因稳定" 是否解耦
# 任务: speed ~ height + drop (两列相关 0.95), bootstrap 200 次
# 对照两种样本量: 全部有落差记录的车 vs n=50 (监督版真实标签集的量级)
import csv, random, math

SRC = 'rcdb_strict.csv'
random.seed(42)

rows = []
with open(SRC) as f:
    for r in csv.DictReader(f):
        try:
            h = float(r['Height']); d = float(r['Drop']); s = float(r['Speed'])
            if h > 0 and d > 0 and s > 0:
                rows.append((h, d, s))
        except (ValueError, TypeError):
            pass
print(f'高度/落差/极速齐全: {len(rows)} 台')

def std(v):
    m = sum(v)/len(v); sd = math.sqrt(sum((x-m)**2 for x in v)/len(v))
    return [(x-m)/sd for x in v], m, sd

def ols2(data):
    # y ~ x1 + x2, 全部标准化后无截距
    x1, _, _ = std([r[0] for r in data])
    x2, _, _ = std([r[1] for r in data])
    y, _, _ = std([r[2] for r in data])
    n = len(data)
    a11 = sum(a*a for a in x1); a12 = sum(a*b for a, b in zip(x1, x2)); a22 = sum(b*b for b in x2)
    b1 = sum(a*c for a, c in zip(x1, y)); b2 = sum(b*c for b, c in zip(x2, y))
    det = a11*a22 - a12*a12
    beta1 = (a22*b1 - a12*b2)/det
    beta2 = (a11*b2 - a12*b1)/det
    pred = [beta1*a + beta2*b for a, b in zip(x1, x2)]
    r2 = 1 - sum((p-c)**2 for p, c in zip(pred, y))/n
    return beta1, beta2, r2

def boot(n_sample, B=200):
    b1s, b2s, sums, r2s = [], [], [], []
    for _ in range(B):
        samp = [random.choice(rows) for _ in range(n_sample)]
        try:
            b1, b2, r2 = ols2(samp)
        except ZeroDivisionError:
            continue
        b1s.append(b1); b2s.append(b2); sums.append(b1+b2); r2s.append(r2)
    def rng(v):
        v = sorted(v)
        return v[int(0.025*len(v))], v[int(0.975*len(v))]
    def mean(v): return sum(v)/len(v)
    flips1 = sum(1 for x in b1s if x*mean(b1s) < 0)
    flips2 = sum(1 for x in b2s if x*mean(b2s) < 0)
    print(f'\n== n={n_sample}, bootstrap {len(b1s)} 次 (标准化系数) ==')
    print(f'  高度系数   均值{mean(b1s):+.3f}  95%区间[{rng(b1s)[0]:+.3f}, {rng(b1s)[1]:+.3f}]  符号翻转 {flips1} 次')
    print(f'  落差系数   均值{mean(b2s):+.3f}  95%区间[{rng(b2s)[0]:+.3f}, {rng(b2s)[1]:+.3f}]  符号翻转 {flips2} 次')
    print(f'  两系数之和 均值{mean(sums):+.3f}  95%区间[{rng(sums)[0]:+.3f}, {rng(sums)[1]:+.3f}]')
    print(f'  R²         均值{mean(r2s):.3f}   95%区间[{rng(r2s)[0]:.3f}, {rng(r2s)[1]:.3f}]')

# 单变量基线
b1, b2, r2_both = ols2(rows)
x1s, _, _ = std([r[0] for r in rows]); ys, _, _ = std([r[2] for r in rows])
r_h = sum(a*b for a, b in zip(x1s, ys))/len(rows)
x2s, _, _ = std([r[1] for r in rows])
r_d = sum(a*b for a, b in zip(x2s, ys))/len(rows)
print(f'只用高度 R²={r_h*r_h:.3f}  只用落差 R²={r_d*r_d:.3f}  两个都用 R²={r2_both:.3f}')
print(f'全样本点估计: 高度系数{b1:+.3f}  落差系数{b2:+.3f}  (标准化)')

boot(len(rows))
boot(50)