因果推断与 A/B 实验：从潜在结果到双重差分

import pandas as pd

df = pd.DataFrame({"group": ["A", "B", "A", "B"],
                   "severity": ["轻", "轻", "重", "重"],
                   "n": [100, 900, 900, 100],
                   "cured": [90, 810, 450, 40]})
df["rate"] = df["cured"] / df["n"]
print(df)
print("合并 A 治愈率:", 540 / 1000)   # 0.54
print("合并 B 治愈率:", 850 / 1000)   # 0.85

P(Y | do(X=1))   vs   P(Y | X=1)
前者：干预 X 后的结果分布（因果）
后者：观察到 X=1 的样本的结果分布（相关）

1. 无干扰（SUTVA）：一个人的结果不受他人干预状态影响
2. 可忽略性（无未观测混杂）：给定协变量 X 后，T 与潜在结果独立
3. 重叠（正值性）：每种协变量组合下都有一定概率被干预

X 干预  →  Y 结果
   ↑          ↑
   └── Z ─────┘      Z 是混杂因子，必须控制
X  →  M  →  Y        M 是中介，控制它会屏蔽部分真实效应
X  →  C  ←  Y        C 是碰撞点，控制它会引入虚假关联

from pgmpy.models import BayesianNetwork
from pgmpy.inference import CausalInference

dag = BayesianNetwork([("Z", "X"), ("Z", "Y"), ("X", "Y")])
print(CausalInference(dag).get_all_backdoor_adjustment_sets("X", "Y"))
# frozenset({frozenset({'Z'})}) —— 控制了 Z 就能识别因果效应

import numpy as np
from sklearn.ensemble import GradientBoostingClassifier

def estimate_ps(X, T):
    """X: 协变量矩阵, T: 干预指示（0/1）"""
    m = GradientBoostingClassifier(n_estimators=200, max_depth=3, random_state=42)
    m.fit(X, T)
    return m.predict_proba(X)[:, 1]

ps = estimate_ps(X, T)
print("倾向得分分布:", np.percentile(ps, [5, 50, 95]).round(3))

from sklearn.neighbors import NearestNeighbors

def match_1to1(ps, T, caliper=0.05):
    treated = np.where(T == 1)[0]
    control = np.where(T == 0)[0]
    nn = NearestNeighbors(n_neighbors=1).fit(ps[control].reshape(-1, 1))
    dist, idx = nn.kneighbors(ps[treated].reshape(-1, 1))
    keep = dist.ravel() <= caliper          # 卡钳：丢弃差距过大的对
    return treated[keep], control[idx.ravel()[keep]]

t_idx, c_idx = match_1to1(ps, T)
print("ATT 估计:", round(Y[t_idx].mean() - Y[c_idx].mean(), 4))

def ipw_ate(Y, T, ps, clip=0.01):
    ps = np.clip(ps, clip, 1 - clip)        # 防极端权重
    w1, w0 = T / ps, (1 - T) / (1 - ps)
    return (w1 * Y).sum() / w1.sum() - (w0 * Y).sum() / w0.sum()

def smd(x, T):                              # 标准化均值差，< 0.1 视为平衡
    m1, m0 = x[T == 1].mean(), x[T == 0].mean()
    s1, s0 = x[T == 1].var(), x[T == 0].var()
    return abs(m1 - m0) / np.sqrt((s1 + s0) / 2)

DID = (Y_treat,post - Y_treat,pre) - (Y_ctrl,post - Y_ctrl,pre)

def did(y_t_post, y_t_pre, y_c_post, y_c_pre):
    return (y_t_post - y_t_pre) - (y_c_post - y_c_pre)

print("DID 效应:", did(120, 100, 115, 105))   # (20) - (10) = 10

import statsmodels.formula.api as smf
res = smf.ols("y ~ treat + post + treat:post", data=df).fit()
print(res.params["treat:post"], res.pvalues["treat:post"])

from scipy.optimize import minimize

def synth_weights(Y_pre_treat, Y_pre_donors):
    """Y_pre_treat: (T,) 处理单位干预前序列; Y_pre_donors: (T, J)"""
    J = Y_pre_donors.shape[1]
    def loss(w):
        w = np.abs(w); w = w / w.sum()
        return np.mean((Y_pre_donors @ w - Y_pre_treat) ** 2)
    res = minimize(loss, np.ones(J) / J, method="Nelder-Mead")
    w = np.abs(res.x)
    return w / w.sum()

1. 明确假设与指标：主指标（唯一）、护栏指标、诊断指标
2. 算样本量：给定基线、最小可检测效应 MDE、alpha、power
3. 随机化与分流：按用户/会话/设备维度，保证一致性
4. 定运行时长：至少覆盖一个完整业务周期（通常 1~2 周）

from scipy.stats import norm

def sample_size_per_group(baseline, mde, alpha=0.05, power=0.8):
    """baseline: 基线转化率; mde: 相对提升（0.05 表示提升 5%）"""
    p1, p2 = baseline, baseline * (1 + mde)
    p_bar = (p1 + p2) / 2
    z_a, z_b = norm.ppf(1 - alpha / 2), norm.ppf(power)
    return int(np.ceil(2 * p_bar * (1 - p_bar) * (z_a + z_b) ** 2 / (p2 - p1) ** 2))

print(sample_size_per_group(0.10, 0.05))   # 每组约 6 万
print(sample_size_per_group(0.10, 0.20))   # 提升到 20%，只需约 4 千

import hashlib

def assign_bucket(user_id, salt, num_buckets=10000):
    """对用户做确定性哈希，保证同一用户始终进同一桶"""
    key = f"{salt}:{user_id}".encode()
    return int(hashlib.md5(key).hexdigest()[:8], 16) % num_buckets

def in_treatment(user_id, salt, ratio=0.5):
    return assign_bucket(user_id, salt) < int(ratio * 10000)

from scipy import stats

def ab_test(conv_a, n_a, conv_b, n_b):
    p_a, p_b = conv_a / n_a, conv_b / n_b
    p_pool = (conv_a + conv_b) / (n_a + n_b)
    se = np.sqrt(p_pool * (1 - p_pool) * (1 / n_a + 1 / n_b))
    z = (p_b - p_a) / se
    return {"lift": (p_b - p_a) / p_a, "z": z,
            "p": 2 * (1 - stats.norm.cdf(abs(z))),
            "ci": ((p_b - p_a) - 1.96 * se, (p_b - p_a) + 1.96 * se)}

print(ab_test(1000, 10000, 1100, 10000))

def bonferroni(p_values, alpha=0.05):
    m = len(p_values)
    return [min(p * m, 1.0) for p in p_values], alpha / m

能随机化吗？
  ├── 能 → A/B 实验（黄金标准）
  └── 不能
        ├── 有前后两期 + 对照组 → DID
        ├── 只有单一大单位 → 合成控制
        └── 横截面 + 丰富协变量 → 倾向得分 / IPW