数据稀疏场景下的因果识别：小样本问题的应对策略

Ⅰ. 小样本因果识别的核心挑战

1.1 统计效能的维度灾难

小样本场景下，因果估计面临三重困境：

以医疗场景为例，某新型抗癌药物III期临床试验中，因患者招募困难，最终仅获得n=68的样本量。此时标准PSM方法无法找到足够匹配的对照样本，ATE估计的方差达到0.42，置信区间宽度超过治疗效果本身的3倍。

1.2 识别假设的脆弱性

因果识别的三大核心假设在小样本下变得异常脆弱：

mermaid流程图：

1.3 实例分析：教育干预项目评估

某非营利组织在贫困山区开展数字化教育干预，仅有15所学校参与实验（处理组5所，对照组10所）。关键挑战包括：

• 基线不平衡：处理组学校平均生源质量比对照组高0.8个标准差
• 流失率高：学期末有效样本仅剩47名学生
• 异质性强：学校间资源配置差异系数达0.65

采用传统线性回归得到ATE估计为0.31(p=0.18)，但安慰剂检验显示有37%的概率获得同等或更大的效应值，结果几乎不可信。

Ⅱ. 应对策略一：贝叶斯因果推断框架

2.1 贝叶斯方法的核心优势

贝叶斯框架通过先验知识引入和后验分布量化不确定性，天然适合小样本场景。其优势体现在：

2.2 层次化贝叶斯模型实现

下面我们实现一个适用于小样本的贝叶斯因果推断模型，包含个体异质性和先验知识融合。

import pymc as pm
import arviz as az
import numpy as np
import pandas as pd
from scipy import stats

class BayesianCausalInference:
    """
    小样本因果推断的贝叶斯实现
    支持层次化建模和先验知识注入
    """
    
    def __init__(self, data, outcome_col, treatment_col, covariate_cols):
        """
        初始化贝叶斯因果推断模型
        
        参数:
        ---------
        data : pd.DataFrame
            输入数据框
        outcome_col : str
            结果变量列名
        treatment_col : str
            处理变量列名(0/1)
        covariate_cols : list
            协变量列名列表
        """
        self.data = data.copy()
        self.outcome_col = outcome_col
        self.treatment_col = treatment_col
        self.covariate_cols = covariate_cols
        self.model = None
        self.trace = None
        
    def build_hierarchical_model(self, 
                                 prior_effect_mu=0.0,
                                 prior_effect_sigma=1.0,
                                 include_individual_heterogeneity=True):
        """
        构建层次化贝叶斯模型
        
        参数:
        ---------
        prior_effect_mu : float
            处理效应先验均值(领域知识)
        prior_effect_sigma : float
            处理效应先验标准差(不确定性)
        include_individual_heterogeneity : bool
            是否包含个体异质性项
        """
        # 数据准备
        y = self.data[self.outcome_col].values
        T = self.data[self.treatment_col].values
        X = self.data[self.covariate_cols].values
        
        # 标准化协变量
        X = (X - X.mean(axis=0)) / X.std(axis=0)
        
        with pm.Model() as self.model:
            # 超先验：控制整体收缩强度
            tau = pm.HalfCauchy('tau', beta=1)
            sigma = pm.HalfCauchy('sigma', beta=1)
            
            # 协变量系数（层次化先验）
            beta_mu = pm.Normal('beta_mu', mu=0, sigma=1)
            beta_sigma = pm.HalfCauchy('beta_sigma', beta=1)
            beta = pm.Normal('beta', mu=beta_mu, sigma=beta_sigma * tau, 
                           shape=X.shape[1])
            
            # 处理效应（信息性先验）
            ate = pm.Normal('ate', mu=prior_effect_mu, sigma=prior_effect_sigma)
            
            # 个体异质性（可选）
            if include_individual_heterogeneity:
                # 假设个体效应围绕ATE正态分布
                cate_sigma = pm.HalfCauchy('cate_sigma', beta=0.5)
                cate_offset = pm.Normal('cate_offset', mu=0, sigma=1, shape=len(y))
                cate = pm.Deterministic('cate', ate + cate_offset * cate_sigma)
            else:
                cate = ate
            
            # 线性预测
            mu = pm.Deterministic('mu', 
                                 T * cate + 
                                 np.dot(X, beta))
            
            # 似然函数
            obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=y)
            
        return self.model
    
    def fit(self, draws=2000, tune=1000, chains=4, target_accept=0.95):
        """
        模型拟合（NUTS采样）
        
        参数:
        ---------
        draws : int
            后验采样次数
        tune : int
            调试阶段迭代次数
        chains : int
            MCMC链数量
        target_accept : float
            接受率目标值（小样本建议>0.9）
        """
        if self.model is None:
            raise ValueError("请先调用build_hierarchical_model构建模型")
            
        # 小样本需要更高的接受率和更多调参
        with self.model:
            self.trace = pm.sample(
                draws=draws,
                tune=tune,
                chains=chains,
                target_accept=target_accept,
                return_inferencedata=True,
                idata_kwargs={"log_likelihood": True}
            )
        
        # 计算关键统计量
        self._compute_posterior_stats()
        
        return self.trace
    
    def _compute_posterior_stats(self):
        """计算后验统计量"""
        if self.trace is None:
            return
            
        # ATE后验统计
        ate_posterior = self.trace.posterior['ate'].values.flatten()
        self.ate_mean = np.mean(ate_posterior)
        self.ate_hdi = az.hdi(ate_posterior, hdi_prob=0.95)
        self.ate_rope_prob = stats.norm.cdf(0.1, loc=self.ate_mean, scale=np.std(ate_posterior)) - \
                            stats.norm.cdf(-0.1, loc=self.ate_mean, scale=np.std(ate_posterior))
        
        print(f"处理效应后验均值: {self.ate_mean:.4f}")
        print(f"95% HDI区间: [{self.ate_hdi[0]:.4f}, {self.ate_hdi[1]:.4f}]")
        print(f"ROPE概率(效应在±0.1内): {self.ate_rope_prob:.2%}")
    
    def diagnose(self):
        """模型诊断"""
        if self.trace is None:
            raise ValueError("请先调用fit方法")
            
        # R-hat诊断
        rhat = az.rhat(self.trace, var_names=['ate', 'sigma'])
        print("\n模型收敛诊断：(R-hat应接近1.0)")
        print(rhat.to_string())
        
        # 有效样本量
        ess = az.ess(self.trace, var_names=['ate'])
        print(f"\nATE有效样本量: {ess['ate'].values[0]:.0f}")
        
        # 绘制后验分布
        az.plot_posterior(self.trace, var_names=['ate'], hdi_prob=0.95)
        
    def predict_counterfactual(self, new_data):
        """
        反事实预测
        
        参数:
        ---------
        new_data : pd.DataFrame
            新观测数据
            
        返回:
        ---------
        dict : 包含事实与反事实预测
        """
        if self.trace is None:
            raise ValueError("请先调用fit方法")
            
        X_new = new_data[self.covariate_cols].values
        X_new = (X_new - X_new.mean(axis=0)) / X_new.std(axis=0)
        
        # 提取后验样本
        beta_samples = self.trace.posterior['beta'].values
        ate_samples = self.trace.posterior['ate'].values
        
        # 计算预测
        n_samples = 1000
        idx = np.random.choice(beta_samples.shape[1], n_samples, replace=False)
        
        y0_pred = np.dot(X_new, beta_samples[0, idx, :].mean(axis=0))  # 未处理
        y1_pred = y0_pred + ate_samples[0, idx].mean()  # 已处理
        
        return {
            'y_factual': y0_pred if new_data[self.treatment_col].iloc[0]==0 else y1_pred,
            'y_counterfactual': y1_pred if new_data[self.treatment_col].iloc[0]==0 else y0_pred,
            'ite': y1_pred - y0_pred
        }

# 实例：教育干预项目贝叶斯分析
# 模拟小样本数据
np.random.seed(42)
n = 50  # 小样本
data = pd.DataFrame({
    'student_id': range(n),
    'baseline_score': np.random.normal(60, 10, n),
    'family_income': np.random.lognormal(10, 0.5, n),
    'treatment': np.random.binomial(1, 0.4, n)  # 40%接受干预
})

# 生成结果（真实ATE=0.35）
true_ate = 0.35
data['outcome'] = 0.5 * data['baseline_score'] + \
                  0.1 * np.log(data['family_income']) + \
                  true_ate * data['treatment'] + \
                  np.random.normal(0, 2, n)

# 应用贝叶斯模型
bci = BayesianCausalInference(
    data=data,
    outcome_col='outcome',
    treatment_col='treatment',
    covariate_cols=['baseline_score', 'family_income']
)

# 构建模型（注入领域知识：预期效应0.3±0.2）
bci.build_hierarchical_model(
    prior_effect_mu=0.3,
    prior_effect_sigma=0.2,
    include_individual_heterogeneity=True
)

# 拟合模型
trace = bci.fit(draws=3000, tune=1500, target_accept=0.95)
bci.diagnose()

# 反事实预测示例
new_student = pd.DataFrame({
    'baseline_score': [65],
    'family_income': [25000],
    'treatment': [0]
})
cf_result = bci.predict_counterfactual(new_student)
print(f"\n新学生反事实预测：")
print(f"不接受干预: {cf_result['y_factual']:.2f}")
print(f"接受干预: {cf_result['y_counterfactual']:.2f}")
print(f"个体处理效应: {cf_result['ite']:.2f}")

2.3 先验分布的校准方法

小样本下先验选择至关重要，我们提供三阶段校准流程：

专家先验量化示例：

def elicit_expert_prior(expert_quantiles):
    """
    将专家分位数意见转化为正态先验参数
    
    参数:
    ---------
    expert_quantiles : dict
        专家提供的分位数，如{'lower':0.1, 'median':0.3, 'upper':0.5}
    """
    # 解决矩匹配问题
    from scipy.optimize import minimize
    
    target_median = expert_quantiles['median']
    target_lower = expert_quantiles['lower']
    target_upper = expert_quantiles['upper']
    
    def objective(params):
        mu, sigma = params
        pred_median = stats.norm.ppf(0.5, loc=mu, scale=sigma)
        pred_lower = stats.norm.ppf(0.25, loc=mu, scale=sigma)
        pred_upper = stats.norm.ppf(0.75, loc=mu, scale=sigma)
        
        return (pred_median - target_median)**2 + \
               (pred_lower - target_lower)**2 + \
               (pred_upper - target_upper)**2
    
    result = minimize(objective, x0=[0.3, 0.1], bounds=[(None, None), (0.01, None)])
    return {'mu': result.x[0], 'sigma': result.x[1]}

# 实际应用
expert_opinion = {'lower': 0.15, 'median': 0.3, 'upper': 0.45}
prior_params = elicit_expert_prior(expert_opinion)
print(f"校准后的先验参数: μ={prior_params['mu']:.3f}, σ={prior_params['sigma']:.3f}")

Ⅲ. 应对策略二：元学习与因果推断融合

3.1 元学习在因果推断中的定位

元学习（Meta-Learning）通过在多个相关任务上学习，使模型能够快速适应新任务。在小样本因果推断中，我们将每个历史实验视为一个"任务"，学习因果效应的先验分布，然后在新的小样本实验中快速调整。

mermaid流程图：

3.2 模型无关元学习（MAML）在CATE估计中的应用

我们实现一个基于MAML的**条件平均处理效应（CATE）**估计框架：

import torch
import torch.nn as nn
from torch.utils.data import Dataset, DataLoader
import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler

class CATEEstimator(nn.Module):
    """CATE估计网络"""
    def __init__(self, input_dim, hidden_dims=[128, 64, 32]):
        super().__init__()
        layers = []
        prev_dim = input_dim
        
        for hidden_dim in hidden_dims:
            layers.extend([
                nn.Linear(prev_dim, hidden_dim),
                nn.ReLU(),
                nn.Dropout(0.3)
            ])
            prev_dim = hidden_dim
        
        # 输出层：预测CATE
        layers.append(nn.Linear(prev_dim, 1))
        self.network = nn.Sequential(*layers)
        
    def forward(self, x):
        return self.network(x)

class MAMLforCausalInference:
    """
    用于因果推断的MAML框架
    支持小样本场景下的快速适应
    """
    
    def __init__(self, input_dim, meta_lr=0.001, inner_lr=0.01):
        self.meta_model = CATEEstimator(input_dim)
        self.meta_optimizer = torch.optim.Adam(
            self.meta_model.parameters(), lr=meta_lr
        )
        self.inner_lr = inner_lr
        self.loss_fn = nn.MSELoss()
        
    def inner_loop_update(self, model, support_set, num_steps=5):
        """
        内循环：在支持集上进行快速梯度更新
        
        参数:
        ---------
        model : nn.Module
            当前任务模型
        support_set : tuple
            (X_support, y_support) 支持集数据
        num_steps : int
            内循环步数
            
        返回:
        ---------
        updated_model : 更新后的模型
        """
        X_support, y_support = support_set
        adapted_model = type(model)(model.network[0].in_features)
        adapted_model.load_state_dict(model.state_dict())
        
        # 内循环优化器
        inner_optimizer = torch.optim.SGD(
            adapted_model.parameters(), lr=self.inner_lr
        )
        
        for step in range(num_steps):
            pred = adapted_model(X_support).squeeze()
            loss = self.loss_fn(pred, y_support)
            inner_optimizer.zero_grad()
            loss.backward()
            inner_optimizer.step()
            
        return adapted_model
    
    def meta_train_step(self, tasks_batch):
        """
        元训练步骤：在任务批次上更新元模型
        
        参数:
        ---------
        tasks_batch : list
            任务列表，每个任务包含支持集和查询集
        """
        meta_loss = 0
        task_losses = []
        
        for task in tasks_batch:
            # 提取支持集和查询集
            support_set = (task['X_support'], task['y_support'])
            query_set = (task['X_query'], task['y_query'])
            
            # 内循环适应
            adapted_model = self.inner_loop_update(
                self.meta_model, support_set, num_steps=5
            )
            
            # 在查询集上评估
            pred_query = adapted_model(query_set[0]).squeeze()
            task_loss = self.loss_fn(pred_query, query_set[1])
            meta_loss += task_loss
            task_losses.append(task_loss.item())
        
        # 元更新
        self.meta_optimizer.zero_grad()
        meta_loss.backward()
        self.meta_optimizer.step()
        
        return np.mean(task_losses)
    
    def adapt_to_new_task(self, new_data, num_steps=10):
        """
        适应新的小样本任务
        
        参数:
        ---------
        new_data : pd.DataFrame
            新实验数据
        num_steps : int
            适应步数
            
        返回:
        ---------
        adapted_model : 适应后的模型
        """
        # 准备数据
        feature_cols = [col for col in new_data.columns 
                       if col not in ['outcome', 'treatment', 'true_cate']]
        X = torch.FloatTensor(new_data[feature_cols].values)
        y = torch.FloatTensor(new_data['true_cate'].values)
        
        # 在新数据上微调
        adapted_model = self.inner_loop_update(
            self.meta_model, (X, y), num_steps=num_steps
        )
        
        return adapted_model
    
    def estimate_cate(self, model, data):
        """
        使用适应后的模型估计CATE
        
        参数:
        ---------
        model : nn.Module
            适应后的模型
        data : pd.DataFrame
            需要预测的数据
            
        返回:
        ---------
        cate_estimates : np.array
            CATE估计值
        """
        feature_cols = [col for col in data.columns 
                       if col not in ['outcome', 'treatment', 'true_cate']]
        X = torch.FloatTensor(data[feature_cols].values)
        model.eval()
        with torch.no_grad():
            cate_preds = model(X).numpy().flatten()
        return cate_preds

# 实例：多任务元学习训练
def simulate_meta_tasks(num_tasks=20, samples_per_task=200, test_samples=50):
    """
    模拟多个因果推断任务（历史实验）
    """
    tasks = []
    
    for task_id in range(num_tasks):
        # 每个任务有不同的真实CATE生成机制
        task_theta = np.random.normal(0, 0.5)
        
        # 生成协变量
        X = np.random.randn(samples_per_task, 10)
        # 生成真实CATE（任务特定）
        true_cate = task_theta + 0.3 * X[:, 0] + 0.2 * X[:, 1]**2
        
        # 生成结果
        treatment = np.random.binomial(1, 0.5, samples_per_task)
        outcome = 0.5 * X[:, 2] + true_cate * treatment + np.random.randn(samples_per_task) * 0.1
        
        # 拆分为支持集和查询集
        split_idx = test_samples
        task_data = {
            'X_support': torch.FloatTensor(X[split_idx:]),
            'y_support': torch.FloatTensor(true_cate[split_idx:]),
            'X_query': torch.FloatTensor(X[:split_idx]),
            'y_query': torch.FloatTensor(true_cate[:split_idx]),
            'task_theta': task_theta
        }
        tasks.append(task_data)
    
    return tasks

# 元训练
maml = MAMLforCausalInference(input_dim=10, meta_lr=0.001, inner_lr=0.01)
tasks = simulate_meta_tasks(num_tasks=50, samples_per_task=100)

# 模拟小样本新任务
new_task_X = np.random.randn(30, 10)  # 仅30个样本！
new_task_theta = np.random.normal(0, 0.5)
new_task_cate = new_task_theta + 0.3 * new_task_X[:, 0] + 0.2 * new_task_X[:, 1]**2
new_task_data = pd.DataFrame(new_task_X, columns=[f'feature_{i}' for i in range(10)])
new_task_data['true_cate'] = new_task_cate

print("开始元训练...")
for epoch in range(100):
    # 随机采样任务批次
    batch_size = 8
    batch_tasks = np.random.choice(tasks, batch_size, replace=False).tolist()
    loss = maml.meta_train_step(batch_tasks)
    
    if epoch % 20 == 0:
        print(f"Epoch {epoch}, Meta Loss: {loss:.4f}")

# 适应新任务
print("\n适应新的小样本任务...")
adapted_model = maml.adapt_to_new_task(new_task_data, num_steps=15)

# 估计CATE
cate_estimates = maml.estimate_cate(adapted_model, new_task_data)

# 评估
true_cates = new_task_data['true_cate'].values
mse = np.mean((cate_estimates - true_cates) ** 2)
correlation = np.corrcoef(cate_estimates, true_cates)[0, 1]

print(f"\n新任务性能:")
print(f"MSE: {mse:.4f}")
print(f"与真实CATE的相关系数: {correlation:.4f}")

# 对比：直接训练的模型
direct_model = CATEEstimator(10)
direct_optimizer = torch.optim.Adam(direct_model.parameters(), lr=0.01)

X_direct = torch.FloatTensor(new_task_X)
y_direct = torch.FloatTensor(new_task_cate)

for step in range(100):
    pred = direct_model(X_direct).squeeze()
    loss = nn.MSELoss()(pred, y_direct)
    direct_optimizer.zero_grad()
    loss.backward()
    direct_optimizer.step()

direct_cates = direct_model(X_direct).detach().numpy().flatten()
direct_mse = np.mean((direct_cates - true_cates) ** 2)
direct_correlation = np.corrcoef(direct_cates, true_cates)[0, 1]

print(f"\n直接训练模型性能:")
print(f"MSE: {direct_mse:.4f}")
print(f"与真实CATE的相关系数: {direct_correlation:.4f}")

print(f"\n元学习相对提升:")
print(f"MSE降低: {(direct_mse - mse) / direct_mse:.2%}")
print(f"相关性提升: {(correlation - direct_correlation) / direct_correlation:.2%}")

Ⅳ. 应对策略三：数据增强与合成控制

4.1 合成控制法的现代演进

合成控制法（Synthetic Control）通过构建未处理单元的加权组合来模拟处理组的反事实状态。在小样本下，我们引入贝叶斯模型平均和非线性组合增强鲁棒性。

4.2 贝叶斯非线性合成控制实现

import numpy as np
import pymc as pm
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel

class BayesianSyntheticControl:
    """
    贝叶斯非线性合成控制
    适用于处理组样本极少的场景
    """
    
    def __init__(self, pre_treatment_data, post_treatment_data, 
                 treated_unit_id, donor_pool_ids):
        """
        初始化合成控制模型
        
        参数:
        ---------
        pre_treatment_data : pd.DataFrame
            处理前数据，列为时间，行为单位
        post_treatment_data : pd.DataFrame
            处理后数据
        treated_unit_id : str/int
            处理组单位ID
        donor_pool_ids : list
            潜在对照组单位ID列表
        """
        self.pre_treatment = pre_treatment_data
        self.post_treatment = post_treatment_data
        self.treated_unit = treated_unit_id
        self.donor_pool = donor_pool_ids
        
        # 数据验证
        self._validate_data()
        
    def _validate_data(self):
        """验证数据结构"""
        assert self.treated_unit in self.pre_treatment.index, "处理单位不在数据中"
        assert all(d in self.pre_treatment.index for d in self.donor_pool), "部分对照单位缺失"
        print(f"数据验证通过：{len(self.donor_pool)}个潜在对照单位")
        
    def fit_bayesian_synthetic(self, model_type='linear', **kwargs):
        """
        拟合贝叶斯合成控制模型
        
        参数:
        ---------
        model_type : str
            'linear'或'nonlinear'
        **kwargs : 模型特定参数
        """
        if model_type == 'linear':
            return self._fit_linear_synthetic()
        elif model_type == 'nonlinear':
            return self._fit_gp_synthetic()
        else:
            raise ValueError("model_type必须为'linear'或'nonlinear'")
    
    def _fit_linear_synthetic(self):
        """线性合成控制（贝叶斯版本）"""
        y_treated_pre = self.pre_treatment.loc[self.treated_unit].values
        X_donors_pre = self.pre_treatment.loc[self.donor_pool].values.T
        
        with pm.Model() as self.linear_model:
            # 权重先验：Dirichlet保证权重非负且和为1
            alpha = np.ones(len(self.donor_pool))  # 均匀先验
            weights = pm.Dirichlet('weights', a=alpha)
            
            # 观测模型
            sigma = pm.HalfNormal('sigma', sigma=1)
            
            # 合成控制预测
            mu = pm.math.dot(X_donors_pre, weights)
            
            # 似然函数
            obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=y_treated_pre)
        
        # 采样
        with self.linear_model:
            self.linear_trace = pm.sample(2000, tune=1000, chains=4)
            
        return self.linear_trace
    
    def _fit_gp_synthetic(self, length_scale=1.0, noise_level=0.1):
        """
        高斯过程非线性合成控制
        
        参数:
        ---------
        length_scale : float
            RBF核的长度尺度
        noise_level : float
            观测噪声水平
        """
        # 准备数据
        y_treated_pre = self.pre_treatment.loc[self.treated_unit].values
        X_donors_pre = self.pre_treatment.loc[self.donor_pool].values
        
        # 使用GP学习非线性权重函数
        kernel = RBF(length_scale=length_scale) + WhiteKernel(noise_level=noise_level)
        self.gp_model = GaussianProcessRegressor(kernel=kernel, alpha=0.1)
        
        # 训练数据：对照组特征 -> 处理组结果
        # 这里将每个时间点的对照组向量作为特征
        self.gp_model.fit(X_donors_pre.T, y_treated_pre)
        
        return self.gp_model
    
    def predict_post_treatment(self, model_type='linear'):
        """
        预测处理后的反事实结果
        
        参数:
        ---------
        model_type : str
            使用的模型类型
            
        返回:
        ---------
        counterfactual_pred : array
            反事实预测值
        uncertainty : array
            预测不确定性
        """
        if model_type == 'linear':
            if not hasattr(self, 'linear_trace'):
                raise ValueError("请先调用fit_bayesian_synthetic")
            
            X_donors_post = self.post_treatment.loc[self.donor_pool].values.T
            
            # 使用后验权重进行预测
            weights_posterior = self.linear_trace.posterior['weights'].values
            # 形状: (chains, draws, n_donors)
            
            # 计算反事实预测分布
            counterfactual_samples = np.tensordot(
                weights_posterior, X_donors_post, axes=([-1], [0])
            )
            
            return {
                'mean': np.mean(counterfactual_samples),
                'hdi_95': az.hdi(counterfactual_samples, hdi_prob=0.95),
                'std': np.std(counterfactual_samples),
                'samples': counterfactual_samples
            }
            
        elif model_type == 'nonlinear':
            if not hasattr(self, 'gp_model'):
                raise ValueError("请先调用fit_bayesian_synthetic")
            
            X_donors_post = self.post_treatment.loc[self.donor_pool].values
            
            # GP预测（包含不确定性）
            pred_mean, pred_std = self.gp_model.predict(
                X_donors_post.T, return_std=True
            )
            
            return {
                'mean': pred_mean,
                'std': pred_std,
                'hdi_95': np.array([pred_mean - 1.96*pred_std, 
                                  pred_mean + 1.96*pred_std])
            }
    
    def estimate_treatment_effect(self, model_type='linear'):
        """
        估计处理效应（ATT）
        
        参数:
        ---------
        model_type : str
            模型类型
            
        返回:
        ---------
        effect_summary : dict
            处理效应统计摘要
        """
        # 实际观测结果
        y_treated_post_actual = self.post_treatment.loc[self.treated_unit].values
        
        # 反事实预测
        cf_pred = self.predict_post_treatment(model_type)
        
        # 处理效应
        treatment_effect = y_treated_post_actual - cf_pred['mean']
        
        # 不确定性传播
        if model_type == 'linear':
            # 从后验样本计算效应分布
            y_actual_expanded = np.expand_dims(y_treated_post_actual, 
                                              (0, 1, 2))
            effect_samples = y_actual_expanded - cf_pred['samples']
            
            return {
                'point_estimate': np.mean(effect_samples),
                'hdi_95': az.hdi(effect_samples, hdi_prob=0.95),
                'probability_positive': np.mean(effect_samples > 0)
            }
        else:
            # GP的解析不确定性
            effect_hdi = np.array([
                y_treated_post_actual - (cf_pred['mean'] + 1.96*cf_pred['std']),
                y_treated_post_actual - (cf_pred['mean'] - 1.96*cf_pred['std'])
            ])
            
            return {
                'point_estimate': treatment_effect,
                'hdi_95': effect_hdi,
                'probability_positive': 1 - stats.norm.cdf(
                    0, loc=treatment_effect, scale=cf_pred['std']
                )
            }

# 实例：加州控烟政策效果评估（小样本场景）
# 模拟数据：只观察3年处理后数据，且只有5个对照州
np.random.seed(42)

# 处理前数据（10年）
time_pre = pd.date_range('1980', '1990', freq='Y')
states = ['California', 'Texas', 'NewYork', 'Florida', 'Illinois', 'Pennsylvania']
treated_state = 'California'

# 生成合成控制数据
pre_data = pd.DataFrame(
    np.random.randn(len(states), len(time_pre)).cumsum(axis=1) + 20,
    index=states,
    columns=time_pre
)
# 模拟加州的真实趋势（略高于合成控制）
pre_data.loc['California'] += np.linspace(0, 0.5, len(time_pre))

# 处理后数据（3年：政策实施）
time_post = pd.date_range('1990', '1993', freq='Y')
post_data = pd.DataFrame(
    np.random.randn(len(states), len(time_post)).cumsum(axis=1) + 25,
    index=states,
    columns=time_post
)
# 政策效应：加州的真实结果下降（政策有效）
post_data.loc['California'] -= np.array([1.5, 1.8, 2.0])

# 应用模型
bsc = BayesianSyntheticControl(
    pre_treatment_data=pre_data,
    post_treatment_data=post_data,
    treated_unit_id=treated_state,
    donor_pool_ids=[s for s in states if s != treated_state]
)

print("=== 线性贝叶斯合成控制 ===")
linear_trace = bsc.fit_bayesian_synthetic(model_type='linear')
linear_effect = bsc.estimate_treatment_effect(model_type='linear')
print(f"ATT点估计: {linear_effect['point_estimate']:.4f}")
print(f"95% HDI: [{linear_effect['hdi_95'][0]:.4f}, {linear_effect['hdi_95'][1]:.4f}]")
print(f"效应为正的概率: {linear_effect['probability_positive']:.2%}")

print("\n=== 高斯过程合成控制 ===")
gp_model = bsc.fit_bayesian_synthetic(model_type='nonlinear')
gp_effect = bsc.estimate_treatment_effect(model_type='nonlinear')
print(f"ATT点估计: {gp_effect['point_estimate']:.4f}")
print(f"95% 置信区间: [{gp_effect['hdi_95'][0]:.4f}, {gp_effect['hdi_95'][1]:.4f}]")
print(f"效应为正的概率: {gp_effect['probability_positive']:.2%}")

# 可视化
import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# 线性模型结果
time_combined = time_pre.union(time_post)
cf_linear = bsc.predict_post_treatment('linear')
ax1.plot(time_pre, pre_data.loc['California'], 'ko-', label='California (Actual)')
ax1.plot(time_post, post_data.loc['California'], 'ro-', label='California (Post)')
ax1.axvline(time_post[0], color='gray', linestyle='--', alpha=0.5)
ax1.fill_between(time_post, 
                 cf_linear['hdi_95'][0], 
                 cf_linear['hdi_95'][1],
                 alpha=0.3, color='blue', label='Synthetic Control 95% HDI')
ax1.set_title('Linear Bayesian Synthetic Control')
ax1.legend()
ax1.set_ylabel('Outcome')

# GP模型结果
cf_gp = bsc.predict_post_treatment('nonlinear')
ax2.plot(time_pre, pre_data.loc['California'], 'ko-')
ax2.plot(time_post, post_data.loc['California'], 'ro-')
ax2.axvline(time_post[0], color='gray', linestyle='--', alpha=0.5)
ax2.fill_between(time_post, 
                 cf_gp['hdi_95'][0], 
                 cf_gp['hdi_95'][1],
                 alpha=0.3, color='green', label='GP Synthetic Control 95% CI')
ax2.set_title('Gaussian Process Synthetic Control')
ax2.legend()

plt.tight_layout()
plt.savefig('synthetic_control_comparison.png', dpi=300)
plt.show()

Ⅴ. 应对策略四：迁移学习与领域适应

5.1 跨域因果知识迁移框架

当目标域数据极度稀疏时，可从源域迁移因果知识。核心挑战在于处理协变量分布偏移和因果机制差异。

mermaid流程图：

5.2 实例：从大规模A/B测试迁移到冷启动业务

场景：电商平台有成熟业务的万级A/B测试数据（源域），需要评估新业务线（目标域）的干预效果，但新业务仅积累200个样本。

5.2.1 领域自适应网络实现

import torch
import torch.nn as nn
import torch.optim as optim
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
import numpy as np
import pandas as pd

class DomainAdversarialNetwork(nn.Module):
    """
    领域对抗网络用于因果推断
    学习领域不变的表示
    """
    
    def __init__(self, input_dim, hidden_dim=64, dropout=0.3):
        super().__init__()
        
        # 特征提取器（共享）
        self.feature_extractor = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.ReLU(),
            nn.Dropout(dropout),
            nn.Linear(hidden_dim, hidden_dim),
            nn.ReLU(),
            nn.Dropout(dropout)
        )
        
        # 预测头（因果效应）
        self.causal_head = nn.Sequential(
            nn.Linear(hidden_dim, hidden_dim // 2),
            nn.ReLU(),
            nn.Linear(hidden_dim // 2, 1)
        )
        
        # 领域分类器（对抗训练）
        self.domain_classifier = nn.Sequential(
            nn.Linear(hidden_dim, hidden_dim // 2),
            nn.ReLU(),
            nn.Linear(hidden_dim // 2, 1),
            nn.Sigmoid()
        )
        
    def forward(self, x, alpha=1.0):
        """
        前向传播
        
        参数:
        ---------
        x : tensor
            输入特征
        alpha : float
            梯度反转层系数
            
        返回:
        ---------
        cate_pred : tensor
            CATE预测
        domain_pred : tensor
            领域预测
        features : tensor
            学习到的特征
        """
        features = self.feature_extractor(x)
        
        # 因果预测
        cate_pred = self.causal_head(features)
        
        # 领域预测（带梯度反转）
        if self.training:
            # 梯度反转操作（手动实现）
            reversed_features = ReverseLayerF.apply(features, alpha)
            domain_pred = self.domain_classifier(reversed_features)
        else:
            domain_pred = self.domain_classifier(features)
            
        return cate_pred, domain_pred, features

class ReverseLayerF(torch.autograd.Function):
    """梯度反转层"""
    @staticmethod
    def forward(ctx, x, alpha):
        ctx.alpha = alpha
        return x.view_as(x)
    
    @staticmethod
    def backward(ctx, grad_output):
        output = grad_output.neg() * ctx.alpha
        return output, None

class TransferCausalLearner:
    """
    迁移学习因果推断主类
    支持源域到目标域的知识迁移
    """
    
    def __init__(self, source_data, target_data, feature_cols,
                 treatment_col, outcome_col):
        """
        初始化迁移学习器
        
        参数:
        ---------
        source_data : pd.DataFrame
            源域数据（大样本）
        target_data : pd.DataFrame
            目标域数据（小样本）
        feature_cols : list
            特征列名
        treatment_col : str
            处理变量列名
        outcome_col : str
            结果变量列名
        """
        self.source_data = source_data
        self.target_data = target_data
        self.feature_cols = feature_cols
        self.treatment_col = treatment_col
        self.outcome_col = outcome_col
        
        # 设备
        self.device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
        
        # 标准化
        self.scaler = StandardScaler()
        all_features = pd.concat([source_data[feature_cols], 
                                 target_data[feature_cols]])
        self.scaler.fit(all_features)
        
    def prepare_data(self, batch_size=64):
        """
        准备数据加载器
        
        参数:
        ---------
        batch_size : int
            批次大小
        """
        # 源域数据
        source_X = self.scaler.transform(self.source_data[self.feature_cols])
        source_T = self.source_data[self.treatment_col].values
        source_Y = self.source_data[self.outcome_col].values
        
        # 目标域数据
        target_X = self.scaler.transform(self.target_data[self.feature_cols])
        target_T = self.target_data[self.treatment_col].values
        target_Y = self.target_data[self.outcome_col].values
        
        # 构建CATE目标（源域用真实差异，目标域用代理）
        source_cate = self._estimate_proxy_cate(self.source_data)
        target_cate = self._estimate_proxy_cate(self.target_data)
        
        # 数据加载器
        source_loader = DataLoader(
            torch.utils.data.TensorDataset(
                torch.FloatTensor(source_X),
                torch.FloatTensor(source_cate),
                torch.zeros(len(source_X))  # 源域标签为0
            ),
            batch_size=batch_size, shuffle=True
        )
        
        target_loader = DataLoader(
            torch.utils.data.TensorDataset(
                torch.FloatTensor(target_X),
                torch.FloatTensor(target_cate),
                torch.ones(len(target_X))   # 目标域标签为1
            ),
            batch_size=batch_size, shuffle=True
        )
        
        return source_loader, target_loader
    
    def _estimate_proxy_cate(self, data):
        """
        估计代理CATE（目标变量对处理的回归系数）
        小样本下使用正则化逻辑回归
        """
        X = data[self.feature_cols]
        T = data[self.treatment_col]
        Y = data[self.outcome_col]
        
        # 简单代理：线性回归系数作为CATE近似
        from sklearn.linear_model import Ridge
        model = Ridge(alpha=1.0)
        model.fit(X, Y)
        
        # 计算预测差异作为代理CATE
        X_treated = X.copy()
        X_treated[self.treatment_col] = 1
        X_control = X.copy()
        X_control[self.treatment_col] = 0
        
        cate_proxy = model.predict(X_treated) - model.predict(X_control)
        return cate_proxy
    
    def train(self, epochs=50, causal_weight=1.0, domain_weight=0.5):
        """
        训练领域对抗网络
        
        参数:
        ---------
        epochs : int
            训练轮数
        causal_weight : float
            因果损失权重
        domain_weight : float
            领域分类损失权重
        """
        input_dim = len(self.feature_cols) + 1  # +1 for treatment
        self.model = DomainAdversarialNetwork(input_dim).to(self.device)
        
        source_loader, target_loader = self.prepare_data()
        
        # 优化器
        optimizer = optim.Adam(self.model.parameters(), lr=0.001)
        
        # 损失函数
        causal_criterion = nn.MSELoss()
        domain_criterion = nn.BCELoss()
        
        self.train_history = []
        
        for epoch in range(epochs):
            self.model.train()
            total_loss = 0
            n_batches = 0
            
            # 迭代源域和目标域
            for (source_batch, target_batch) in zip(source_loader, target_loader):
                # 源域数据
                source_x, source_cate, source_domain = source_batch
                source_x = source_x.to(self.device)
                source_cate = source_cate.to(self.device).unsqueeze(1)
                source_domain = source_domain.to(self.device).unsqueeze(1)
                
                # 目标域数据
                target_x, target_cate, target_domain = target_batch
                target_x = target_x.to(self.device)
                target_cate = target_cate.to(self.device).unsqueeze(1)
                target_domain = target_domain.to(self.device).unsqueeze(1)
                
                # 合并数据
                combined_x = torch.cat([source_x, target_x], dim=0)
                combined_cate = torch.cat([source_cate, target_cate], dim=0)
                combined_domain = torch.cat([source_domain, target_domain], dim=0)
                
                # 前向传播
                cate_pred, domain_pred, _ = self.model(combined_x, alpha=1.0)
                
                # 计算损失
                causal_loss = causal_criterion(cate_pred, combined_cate)
                domain_loss = domain_criterion(domain_pred, combined_domain)
                
                # 总损失
                loss = causal_weight * causal_loss + domain_weight * domain_loss
                
                # 反向传播
                optimizer.zero_grad()
                loss.backward()
                optimizer.step()
                
                total_loss += loss.item()
                n_batches += 1
            
            avg_loss = total_loss / n_batches
            self.train_history.append(avg_loss)
            
            if epoch % 10 == 0:
                print(f"Epoch {epoch}: Loss = {avg_loss:.4f}")
        
        print("训练完成")
    
    def estimate_target_ate(self):
        """
        估计目标域的平均处理效应
        """
        self.model.eval()
        
        # 目标域数据
        target_X = self.scaler.transform(self.target_data[self.feature_cols])
        target_T = self.target_data[self.treatment_col].values
        
        # 构建输入（包含处理变量）
        X_control = np.column_stack([target_X, np.zeros_like(target_T)])
        X_treated = np.column_stack([target_X, np.ones_like(target_T)])
        
        with torch.no_grad():
            cate_control, _, _ = self.model(
                torch.FloatTensor(X_control).to(self.device)
            )
            cate_treated, _, _ = self.model(
                torch.FloatTensor(X_treated).to(self.device)
            )
        
        # ATE估计
        ate = (cate_treated - cate_control).cpu().numpy().mean()
        return ate
    
    def importance_weighting_adjustment(self):
        """
        重要性权重调整（处理协变量偏移）
        返回权重用于后续加权估计
        """
        # 计算密度比 w(x) = p_target(x) / p_source(x)
        source_X = self.scaler.transform(self.source_data[self.feature_cols])
        target_X = self.scaler.transform(self.target_data[self.feature_cols])
        
        # 使用逻辑回归估计密度比
        X_combined = np.vstack([source_X, target_X])
        y_domain = np.hstack([
            np.zeros(len(source_X)), 
            np.ones(len(target_X))
        ])
        
        # 训练领域判别器
        domain_clf = LogisticRegression(max_iter=1000)
        domain_clf.fit(X_combined, y_domain)
        
        # 计算权重: w(x) = p(target) / p(source) = exp(logit) * (n_target/n_source)
        source_pred_logit = domain_clf.decision_function(source_X)
        weights = np.exp(source_pred_logit) * (len(target_X) / len(source_X))
        
        # 截断权重防止方差爆炸
        weights = np.clip(weights, 0.1, 10)
        
        return weights

# 实例：电商营销效果迁移
# 模拟源域（成熟业务）和目标域（新业务）数据
np.random.seed(42)

# 源域数据：10000样本
n_source = 10000
source_data = pd.DataFrame({
    'user_age': np.random.normal(35, 10, n_source),
    'purchase_history': np.random.lognormal(5, 1, n_source),
    'browsing_time': np.random.gamma(2, 30, n_source),
    'treatment': np.random.binomial(1, 0.5, n_source)
})

# 源域真实CATE：复杂非线性
source_cate = 0.3 + \
              0.1 * (source_data['user_age'] - 35) / 10 + \
              0.05 * np.log(source_data['purchase_history']) + \
              0.01 * source_data['browsing_time']
source_data['outcome'] = 0.5 * source_data['purchase_history'] + \
                         source_cate * source_data['treatment'] + \
                         np.random.normal(0, 2, n_source)

# 目标域数据：仅200样本（存在分布偏移）
n_target = 200
target_data = pd.DataFrame({
    'user_age': np.random.normal(28, 8, n_target),  # 更年轻
    'purchase_history': np.random.lognormal(4, 0.8, n_target),  # 消费更低
    'browsing_time': np.random.gamma(1.5, 25, n_target),
    'treatment': np.random.binomial(1, 0.4, n_target)  # 不同处理分配
})

# 目标域CATE（机制相似但参数不同）
target_cate = 0.25 + \
              0.08 * (target_data['user_age'] - 28) / 8 + \
              0.06 * np.log(target_data['purchase_history']) + \
              0.015 * target_data['browsing_time']
target_data['outcome'] = 0.4 * target_data['purchase_history'] + \
                         target_cate * target_data['treatment'] + \
                         np.random.normal(0, 2.5, n_target)

feature_cols = ['user_age', 'purchase_history', 'browsing_time']

# 迁移学习
transfer_learner = TransferCausalLearner(
    source_data=source_data,
    target_data=target_data,
    feature_cols=feature_cols,
    treatment_col='treatment',
    outcome_col='outcome'
)

# 训练
print("开始领域对抗训练...")
transfer_learner.train(epochs=30, causal_weight=1.0, domain_weight=0.5)

# 估计目标域ATE
target_ate = transfer_learner.estimate_target_ate()
print(f"\n迁移学习估计的目标域ATE: {target_ate:.4f}")

# 基线：仅在目标域训练
baseline_model = DomainAdversarialNetwork(len(feature_cols) + 1)
baseline_criterion = nn.MSELoss()
baseline_optimizer = optim.Adam(baseline_model.parameters(), lr=0.001)

target_X_scaled = transfer_learner.scaler.transform(target_data[feature_cols])
target_cate_proxy = transfer_learner._estimate_proxy_cate(target_data)

baseline_X = torch.FloatTensor(np.column_stack([
    target_X_scaled, 
    target_data['treatment'].values
]))
baseline_y = torch.FloatTensor(target_cate_proxy)

for epoch in range(100):
    baseline_pred, _, _ = baseline_model(baseline_X)
    loss = baseline_criterion(baseline_pred.squeeze(), baseline_y)
    baseline_optimizer.zero_grad()
    loss.backward()
    baseline_optimizer.step()

baseline_ate = (baseline_pred[:, 1] - baseline_pred[:, 0]).mean().item()
print(f"仅目标域训练的ATE: {baseline_ate:.4f}")

# 真实效应对比
true_target_ate = target_cate.mean()
print(f"真实目标域ATE: {true_target_ate:.4f}")
print(f"迁移学习误差: {abs(target_ate - true_target_ate):.4f}")
print(f"基线方法误差: {abs(baseline_ate - true_target_ate):.4f}")
print(f"迁移学习相对提升: {abs(baseline_ate - true_target_ate) / abs(target_ate - true_target_ate) - 1:.2%}")

# 重要性权重
importance_weights = transfer_learner.importance_weighting_adjustment()
print(f"\n重要性权重统计：")
print(f"均值: {importance_weights.mean():.2f}")
print(f"标准差: {importance_weights.std():.2f}")

Ⅵ. 应对策略五：主动学习与实验设计优化

6.1 信息价值驱动的实验设计

当数据获取成本高昂时，主动学习通过智能选择最有信息量的样本进行干预，最大化统计效能。核心思想是将因果估计的不确定性转化为采样策略。

6.2 贝叶斯主动学习实现

import numpy as np
import pymc as pm
from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.model_selection import train_test_split
import pandas as pd
from scipy.stats import entropy

class BayesianActiveLearning:
    """
    贝叶斯主动学习框架
    用于小样本实验的批次选择
    """
    
    def __init__(self, pool_data, feature_cols, budget=50):
        """
        初始化主动学习器
        
        参数:
        ---------
        pool_data : pd.DataFrame
            未标注样本池
        feature_cols : list
            特征列名
        budget : int
            总预算（最大样本数）
        """
        self.pool_data = pool_data.copy()
        self.feature_cols = feature_cols
        self.budget = budget
        self.labeled_indices = []
        self.model = None
        self.trace = None
        
    def initialize_random(self, n_initial=10):
        """
        随机初始化标注集
        
        参数:
        ---------
        n_initial : int
            初始标注样本数
        """
        initial_idx = np.random.choice(
            len(self.pool_data), n_initial, replace=False
        )
        self.labeled_indices = initial_idx.tolist()
        print(f"随机初始化 {n_initial} 个样本")
        return initial_idx
    
    def build_bayesian_gp_model(self, X, y):
        """
        构建贝叶斯高斯过程模型
        
        参数:
        ---------
        X : array
            标注特征
        y : array
            标注结果（处理效应）
        """
        with pm.Model() as self.model:
            # 核函数
            ℓ = pm.Gamma('ℓ', alpha=2, beta=1)  # 长度尺度
            η = pm.HalfCauchy('η', beta=5)      # 信号方差
            
            cov = η**2 * pm.gp.cov.Matern32(X.shape[1], ℓ) + \
                  pm.gp.cov.WhiteNoise(0.1)
            
            # 高斯过程
            gp = pm.gp.Marginal(cov_func=cov)
            
            # 先验
            f = gp.marginal_likelihood('f', X=X, y=y, noise=0.1)
        
        return self.model
    
    def fit(self, X, y, draws=1000, tune=500):
        """
        拟合模型
        
        参数:
        ---------
        X : array
            训练特征
        y : array
            训练标签
        draws : int
            后验采样数
        tune : int
            调参迭代数
        """
        self.build_bayesian_gp_model(X, y)
        
        with self.model:
            self.trace = pm.sample(draws=draws, tune=tune, chains=4)
        
        return self.trace
    
    def acquire_points(self, acquisition_type='bald', n_points=5):
        """
        选择下一个批次样本
        
        参数:
        ---------
        acquisition_type : str
            采集函数类型：'bald', 'max_variance', 'thompson'
        n_points : int
            本次选择样本数
        """
        # 未标注池
        unlabeled_idx = np.setdiff1d(
            range(len(self.pool_data)), self.labeled_indices
        )
        
        if len(unlabeled_idx) == 0:
            return []
        
        X_pool = self.pool_data.iloc[unlabeled_idx][self.feature_cols].values
        
        if acquisition_type == 'bald':
            # 贝叶斯主动学习分歧(BALD)
            scores = self._bald_score(X_pool)
        elif acquisition_type == 'max_variance':
            # 最大后验方差
            scores = self._predictive_variance(X_pool)
        elif acquisition_type == 'thompson':
            # Thompson采样
            scores = self._thompson_sampling(X_pool)
        else:
            raise ValueError("不支持的采集函数")
        
        # 选择top-k
        selected_idx = unlabeled_idx[np.argsort(scores)[-n_points:]]
        self.labeled_indices.extend(selected_idx.tolist())
        
        return selected_idx
    
    def _predictive_variance(self, X_pool):
        """
        计算预测方差（不确定性采样）
        """
        if self.trace is None:
            return np.random.random(len(X_pool))
        
        # 从后验采样核参数
        ℓ_samples = self.trace.posterior['ℓ'].values.flatten()
        η_samples = self.trace.posterior['η'].values.flatten()
        
        variances = []
        for ℓ, η in zip(ℓ_samples[::10], η_samples[::10]):  # 子采样加速
            kernel = η**2 * RBF(length_scale=ℓ) + WhiteKernel(0.1)
            gp = GaussianProcessRegressor(kernel=kernel, alpha=0.1)
            
            # 使用已标注数据训练GP
            X_labeled = self.pool_data.iloc[self.labeled_indices][self.feature_cols].values
            y_labeled = self.pool_data.iloc[self.labeled_indices]['observed_effect'].values
            
            gp.fit(X_labeled, y_labeled)
            
            # 预测方差
            _, std = gp.predict(X_pool, return_std=True)
            variances.append(std**2)
        
        return np.mean(variances, axis=0)
    
    def _bald_score(self, X_pool):
        """
        计算BALD分数（贝叶斯主动学习分歧）
        
        BALD = H[y|x,D] - E_{w}[H[y|x,w]]
        """
        if self.trace is None:
            return np.random.random(len(X_pool))
        
        # 预测分布的熵
        predictive_entropy = self._predictive_entropy(X_pool)
        
        # 期望的似然熵（近似）
        expected_likelihood_entropy = self._expected_likelihood_entropy(X_pool)
        
        bald_score = predictive_entropy - expected_likelihood_entropy
        return bald_score
    
    def _predictive_entropy(self, X_pool):
        """计算预测分布的熵"""
        # 使用后验预测采样
        effect_samples = self._sample_effects(X_pool, n_samples=200)
        
        # 离散化计算熵
        entropies = []
        for i in range(X_pool.shape[0]):
            hist, _ = np.histogram(effect_samples[:, i], bins=20, density=True)
            hist = hist[hist > 0]  # 避免log(0)
            entropies.append(-np.sum(hist * np.log(hist)))
        
        return np.array(entropies)
    
    def _sample_effects(self, X_pool, n_samples=100):
        """从后验预测采样效应"""
        ℓ_samples = self.trace.posterior['ℓ'].values.flatten()[:n_samples]
        η_samples = self.trace.posterior['η'].values.flatten()[:n_samples]
        
        all_samples = []
        for ℓ, η in zip(ℓ_samples, η_samples):
            kernel = η**2 * RBF(length_scale=ℓ) + WhiteKernel(0.1)
            gp = GaussianProcessRegressor(kernel=kernel, alpha=0.1)
            
            X_labeled = self.pool_data.iloc[self.labeled_indices][self.feature_cols].values
            y_labeled = self.pool_data.iloc[self.labeled_indices]['observed_effect'].values
            
            gp.fit(X_labeled, y_labeled)
            
            # 采样而非点预测
            samples = gp.sample_y(X_pool, n_samples=1, random_state=None)
            all_samples.append(samples.flatten())
        
        return np.array(all_samples)
    
    def run_active_experiment(self, acquisition_type='bald', batch_size=5):
        """
        运行主动学习实验
        
        参数:
        ---------
        acquisition_type : str
            采集函数类型
        batch_size : int
            每轮选择样本数
        """
        results = []
        
        # 初始化
        self.initialize_random(n_initial=10)
        
        while len(self.labeled_indices) < self.budget:
            # 获取标注数据
            labeled_data = self.pool_data.iloc[self.labeled_indices]
            X_labeled = labeled_data[self.feature_cols].values
            y_labeled = labeled_data['observed_effect'].values
            
            # 拟合模型
            self.fit(X_labeled, y_labeled)
            
            # 评估当前ATE估计
            current_ate = self._estimate_ate()
            ate_error = abs(current_ate - self._true_ate())
            
            results.append({
                'n_samples': len(self.labeled_indices),
                'ate_error': ate_error,
                'acquisition_type': acquisition_type
            })
            
            print(f"已标注 {len(self.labeled_indices)} 样本，ATE误差: {ate_error:.4f}")
            
            # 选择下一批次
            if len(self.labeled_indices) + batch_size <= self.budget:
                new_indices = self.acquire_points(
                    acquisition_type, batch_size
                )
                
                # 模拟标注（在实际中需要真实实验）
                self._simulate_labeling(new_indices)
        
        return pd.DataFrame(results)
    
    def _simulate_labeling(self, indices):
        """模拟标注过程（真实场景中替换为实际实验）"""
        # 在池中标记这些样本为"已观测"
        # 这里假设我们有一个潜在的标签生成机制
        pass
    
    def _estimate_ate(self):
        """基于当前模型估计ATE"""
        if self.trace is None:
            return 0
        
        # 对池中所有样本预测CATE并平均
        X_pool = self.pool_data[self.feature_cols].values
        effect_samples = self._sample_effects(X_pool, n_samples=100)
        return np.mean(effect_samples)
    
    def _true_ate(self):
        """返回真实ATE（用于模拟评估）"""
        return self.pool_data['true_effect'].mean()

# 实例：医疗试验主动学习
# 模拟患者池
np.random.seed(42)
n_pool = 1000
patient_pool = pd.DataFrame({
    'age': np.random.normal(50, 15, n_pool),
    'severity': np.random.beta(2, 5, n_pool),  # 疾病严重程度
    'comorbidity': np.random.poisson(1, n_pool),
    'biomarker': np.random.normal(0, 1, n_pool)
})

# 生成真实个体处理效应（非线性）
true_cate = 0.2 + \
            0.1 * patient_pool['severity'] + \
            0.05 * (patient_pool['age'] - 50) / 15 + \
            0.08 * patient_pool['biomidity'] * patient_pool['biomarker'] + \
            np.random.normal(0, 0.05, n_pool)

patient_pool['true_effect'] = true_cate
patient_pool['observed_effect'] = np.nan  # 初始未观测

# 主动学习实验
budget = 50
acquisition_functions = ['max_variance', 'bald', 'thompson']
results_comparison = []

for acq_fn in acquisition_functions:
    print(f"\n=== 使用 {acq_fn.upper()} 采集函数 ===")
    al = BayesianActiveLearning(patient_pool, 
                                feature_cols=['age', 'severity', 'comorbidity', 'biomarker'],
                                budget=budget)
    
    # 运行实验
    result_df = al.run_active_experiment(acquisition_type=acq_fn, batch_size=5)
    result_df['acq_fn'] = acq_fn
    results_comparison.append(result_df)

# 比较不同策略
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10, 6))
for acq_fn in acquisition_functions:
    df = pd.concat(results_comparison)
    df[df['acq_fn'] == acq_fn].plot(
        x='n_samples', y='ate_error', label=acq_fn, ax=ax
    )

ax.set_xlabel('标注样本数')
ax.set_ylabel('ATE估计误差')
ax.set_title('主动学习采集函数比较')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('active_learning_comparison.png', dpi=300)
plt.show()

# 随机基准
random_errors = []
for _ in range(10):
    al_random = BayesianActiveLearning(patient_pool, 
                                       ['age', 'severity', 'comorbidity', 'biomarker'],
                                       budget=50)
    al_random.initialize_random(n_initial=10)
    while len(al_random.labeled_indices) < 50:
        new_idx = np.random.choice(
            np.setdiff1d(range(n_pool), al_random.labeled_indices),
            5, replace=False
        )
        al_random.labeled_indices.extend(new_idx)
    
    X_labeled = patient_pool.iloc[al_random.labeled_indices][['age', 'severity', 'comorbidity', 'biomarker']].values
    y_labeled = patient_pool.iloc[al_random.labeled_indices]['observed_effect'].values
    al_random.fit(X_labeled, y_labeled)
    random_errors.append(abs(al_random._estimate_ate() - al_random._true_ate()))

print(f"\n随机基准最终误差: {np.mean(random_errors):.4f} ± {np.std(random_errors):.4f}")

# 最优策略最终误差
bald_final = results_comparison[1][results_comparison[1]['n_samples'] == 50]['ate_error'].iloc[0]
print(f"BALD策略最终误差: {bald_final:.4f}")
print(f"相对改进: {(np.mean(random_errors) - bald_final) / np.mean(random_errors):.2%}")

Ⅶ. 综合实战案例：跨境电商用户激励策略评估

7.1 业务背景与挑战

某跨境电商平台在**新兴市场（东南亚）**推出新用户激励计划，面临严峻挑战：

• 数据极度稀疏：上线仅2周，累计用户280人，其中处理组95人
• 高维混杂因素：用户来自5个国家，设备类型、支付方式等30+维度
• 动态环境：汇率波动、政策变化导致时间趋势非平稳
• 成本敏感：每个激励成本$8，需快速验证ROI

目标：准确估计激励对用户首单转化率（二值结果）的因果效应，误差需控制在±2%以内。

7.2 多策略融合解决方案

我们整合前述五种策略，构建层次化融合框架：

mermaid流程图：

7.3 完整代码实现与部署

import pandas as pd
import numpy as np
import pymc as pm
import torch
from sklearn.preprocessing import StandardScaler, LabelEncoder
from sklearn.model_selection import train_test_split
import warnings
warnings.filterwarnings('ignore')

class CrossBorderCausalSystem:
    """
    跨境电商因果推断系统
    集成多策略应对小样本问题
    """
    
    def __init__(self, source_market_data, target_market_data):
        """
        初始化系统
        
        参数:
        ---------
        source_market_data : pd.DataFrame
            成熟市场数据（源域）
        target_market_data : pd.DataFrame
            新兴市场数据（目标域）
        """
        self.source = source_market_data
        self.target = target_market_data
        
        # 特征工程
        self.feature_cols = self._engineer_features()
        
        # 标准化器
        self.scaler = StandardScaler()
        self.scaler.fit(pd.concat([
            self.source[self.feature_cols],
            self.target[self.feature_cols]
        ]))
        
        # 结果存储
        self.models = {}
        self.estimates = {}
        
    def _engineer_features(self):
        """跨市场特征工程"""
        # 对类别变量编码
        le_country = LabelEncoder()
        combined_country = pd.concat([
            self.source['country'], 
            self.target['country']
        ])
        le_country.fit(combined_country)
        
        self.source['country_encoded'] = le_country.transform(self.source['country'])
        self.target['country_encoded'] = le_country.transform(self.target['country'])
        
        # 设备类型编码
        le_device = LabelEncoder()
        combined_device = pd.concat([
            self.source['device_type'], 
            self.target['device_type']
        ])
        le_device.fit(combined_device)
        
        self.source['device_encoded'] = le_device.transform(self.source['device_type'])
        self.transfer.py
        self.target['device_encoded'] = le_device.transform(self.target['device_type'])
        
        return ['country_encoded', 'device_encoded', 'days_since_reg', 
                'browsing_pages', 'cart_items', 'referral_score']
    
    def phase1_data_augmentation(self):
        """
        阶段1：数据增强
        - 基于源域分布的SMOTE变体
        - 时间序列平滑
        """
        from imblearn.over_sampling import SMOTE
        
        # 分离处理组和对照组
        source_treated = self.source[self.source['treatment'] == 1]
        source_control = self.source[self.source['treatment'] == 0]
        
        print(f"源域处理组: {len(source_treated)}, 对照组: {len(source_control)}")
        
        # 对少数类（转化用户）进行过采样
        X_features = self.source[self.feature_cols]
        y_outcome = self.source['converted']
        
        smote = SMOTE(k_neighbors=min(3, len(source_treated)-1), random_state=42)
        X_res, y_res = smote.fit_resample(X_features, y_outcome)
        
        self.augmented_source = pd.DataFrame(X_res, columns=self.feature_cols)
        self.augmented_source['converted'] = y_res
        self.augmented_source['treatment'] = np.random.binomial(1, 0.5, len(self.augmented_source))
        
        print(f"增强后源域样本: {len(self.augmented_source)}")
        return self.augmented_source
    
    def phase2_transfer_learning(self):
        """
        阶段2：迁移学习
        - 领域对抗训练
        - 重要性权重调整
        """
        # 使用增强后的源域
        transfer_data = pd.concat([
            self.augmented_source.assign(domain=0),
            self.target.assign(domain=1)
        ])
        
        # 构建迁移学习模型
        from transfer import TransferCausalLearner
        
        self.transfer_model = TransferCausalLearner(
            source_data=self.augmented_source,
            target_data=self.target,
            feature_cols=self.feature_cols,
            treatment_col='treatment',
            outcome_col='converted'
        )
        
        # 训练
        self.transfer_model.train(epochs=30, causal_weight=1.0, domain_weight=0.6)
        
        # 计算重要性权重
        self.importance_weights = self.transfer_model.importance_weighting_adjustment()
        
        print("迁移学习完成")
        return self.transfer_model
    
    def phase3_bayesian_modeling(self):
        """
        阶段3：贝叶斯层次化建模
        - 融合迁移学习特征
        - 国家层级随机效应
        """
        # 准备数据
        target_X = self.scaler.transform(self.target[self.feature_cols])
        target_T = self.target['treatment'].values
        target_Y = self.target['converted'].values
        
        # 国家层级编码
        country_idx = self.target['country_encoded'].values
        
        with pm.Model() as self.bayesian_model:
            # 超先验
            mu_ate = pm.Normal('mu_ate', mu=0.03, sigma=0.01)  # 预期效应3%
            sigma_country = pm.HalfNormal('sigma_country', sigma=0.02)
            
            # 国家层级效应
            n_countries = len(self.target['country_encoded'].unique())
            country_effect = pm.Normal('country_effect', 
                                       mu=0, 
                                       sigma=sigma_country, 
                                       shape=n_countries)
            
            # 个体效应
            individual_ate = pm.Normal('individual_ate', 
                                       mu=mu_ate + country_effect[country_idx],
                                       sigma=0.015)
            
            # 倾向得分模型（用于调整）
            beta_ps = pm.Normal('beta_ps', mu=0, sigma=1, shape=len(self.feature_cols))
            ps = pm.math.sigmoid(pm.math.dot(target_X, beta_ps))
            
            # 结果模型（考虑迁移特征）
            beta_outcome = pm.Normal('beta_outcome', mu=0, sigma=1, 
                                   shape=len(self.feature_cols))
            
            # 逆概率加权
            ipw = target_T / ps + (1 - target_T) / (1 - ps)
            
            # 似然函数（加权）
            mu_y = pm.math.sigmoid(
                pm.math.dot(target_X, beta_outcome) + individual_ate * target_T
            )
            
            obs = pm.Bernoulli('obs', p=mu_y, observed=target_Y, 
                             size=len(target_Y))  # 重要：传递权重
            
            # 后验
            self.trace = pm.sample(2000, tune=1000, chains=4, target_accept=0.95)
        
        # 提取ATE
        ate_samples = self.trace.posterior['individual_ate'].values.flatten()
        self.estimates['bayes_ate'] = {
            'mean': np.mean(ate_samples),
            'hdi_95': pm.hdi(ate_samples, hdi_prob=0.95),
            'prob_positive': np.mean(ate_samples > 0)
        }
        
        return self.estimates['bayes_ate']
    
    def phase4_synthetic_control_validation(self):
        """
        阶段4：合成控制验证
        - 使用时间序列数据验证趋势
        """
        # 按天聚合数据
        daily_target = self.target.groupby(['date', 'treatment']).agg({
            'converted': 'mean',
            'user_id': 'count'
        }).reset_index()
        
        # 构建合成控制（处理组 vs 对照组）
        treated_daily = daily_target[daily_target['treatment'] == 1].pivot(
            index='date', columns='treatment', values='converted'
        )
        control_daily = daily_target[daily_target['treatment'] == 0].pivot(
            index='date', columns='treatment', values='converted'
        )
        
        # 贝叶斯结构时间序列
        with pm.Model() as self.sts_model:
            # 局部线性趋势
            trend = pm.GaussianRandomWalk('trend', sigma=0.01)
            
            # 季节效应（7天）
            season = pm.Normal('season', mu=0, sigma=0.02, shape=7)
            
            # 合成控制预测
            control_obs = control_daily.values.flatten()
            pred_control = trend + season[:len(control_obs)]
            
            # 观测模型
            obs = pm.Normal('obs', mu=pred_control, sigma=0.02, 
                           observed=control_obs)
            
            self.sts_trace = pm.sample(1000, tune=500)
        
        # 外推预测处理组反事实
        trend_pred = self.sts_trace.posterior['trend'].values[:, -1].mean()
        season_pred = self.sts_trace.posterior['season'].values[:, :7].mean(axis=0)
        
        # 比较实际与预测
        treated_actual = treated_daily.values.flatten()
        treated_counterfactual = trend_pred + season_pred[:len(treated_actual)]
        
        self.estimates['sc_validation'] = {
            'actual_mean': treated_actual.mean(),
            'counterfactual_mean': treated_counterfactual.mean(),
            'effect': treated_actual.mean() - treated_counterfactual.mean()
        }
        
        return self.estimates['sc_validation']
    
    def phase5_active_learning_suggestion(self, pool_data, n_suggest=30):
        """
        阶段5：主动学习建议
        - 识别下一批最信息量用户
        """
        from active_learning import BayesianActiveLearning
        
        # 构建预留池（未处理用户）
        candidate_pool = pool_data[pool_data['treatment'] == 0].copy()
        
        # 添加代理效应（迁移学习预测）
        candidate_pool['proxy_effect'] = self.transfer_model.estimate_cate(
            candidate_pool
        )
        
        al = BayesianActiveLearning(
            data=candidate_pool,
            feature_cols=self.feature_cols,
            budget=n_suggest
        )
        
        # 使用BALD选择
        suggested_idx = al.acquire_points(
            acquisition_type='bald', 
            n_points=n_suggest
        )
        
        self.recommended_users = candidate_pool.iloc[suggested_idx]
        
        print(f"建议下一批干预 {len(suggested_idx)} 个高价值用户")
        return self.recommended_users
    
    def generate_report(self):
        """
        生成综合报告
        """
        report = f"""
        =========================================
        跨境电商用户激励策略因果效应评估报告
        =========================================
        
        数据概况:
        - 目标市场样本: {len(self.target)} (处理组: {self.target['treatment'].sum()})
        - 源市场增强样本: {len(self.augmented_source)}
        - 特征维度: {len(self.feature_cols)}
        
        多策略估计结果:
        
        Ⅰ. 迁移学习ATE: {self.estimates['transfer_ate']:.4f}
        
        Ⅱ. 贝叶斯层次模型:
           - 后验均值: {self.estimates['bayes_ate']['mean']:.4f}
           - 95% HDI: [{self.estimates['bayes_ate']['hdi_95'][0]:.4f}, 
                      {self.estimates['bayes_ate']['hdi_95'][1]:.4f}]
           - 效应为正概率: {self.estimates['bayes_ate']['prob_positive']:.2%}
        
        Ⅲ. 合成控制验证:
           - 实际转化率: {self.estimates['sc_validation']['actual_mean']:.4f}
           - 反事实预测: {self.estimates['sc_validation']['counterfactual_mean']:.4f}
           - 效应一致性: {self.estimates['sc_validation']['effect']:.4f}
        
        Ⅳ. 建议行动:
           - 推荐干预用户数: {len(self.recommended_users)}
           - 预期增量转化: {self.recommended_users['proxy_effect'].sum():.2f}
           - 预计ROI: {self.recommended_users['proxy_effect'].sum() * 50 / (len(self.recommended_users) * 8):.2f}
        
        结论:
        综合多策略估计，激励策略ATE为 {self.estimates['bayes_ate']['mean']:.2%}，
        显著为正的概率 {self.estimates['bayes_ate']['prob_positive']:.0%}，
        建议扩大实验规模至 {len(self.recommended_users)} 用户。
        """
        
        print(report)
        return report

# 数据准备
np.random.seed(42)

# 模拟源域数据（成熟市场：美国）
n_source = 15000
source_data = pd.DataFrame({
    'country': ['USA'] * n_source,
    'device_type': np.random.choice(['iOS', 'Android', 'Web'], n_source, p=[0.4, 0.35, 0.25]),
    'days_since_reg': np.random.exponential(30, n_source),
    'browsing_pages': np.random.negative_binomial(5, 0.3, n_source),
    'cart_items': np.random.poisson(2, n_source),
    'referral_score': np.random.beta(2, 8, n_source),
    'treatment': np.random.binomial(1, 0.5, n_source),
    'converted': np.random.binomial(1, 0.15, n_source)
})

# 模拟目标域数据（新兴市场：东南亚）
countries = ['ID', 'TH', 'VN', 'MY', 'PH']
n_target = 280
target_data = pd.DataFrame({
    'country': np.random.choice(countries, n_target, p=[0.3, 0.25, 0.2, 0.15, 0.1]),
    'device_type': np.random.choice(['Android', 'Web'], n_target, p=[0.7, 0.3]),  # Android主导
    'days_since_reg': np.random.exponential(5, n_target),  # 新用户更多
    'browsing_pages': np.random.negative_binomial(3, 0.4, n_target),
    'cart_items': np.random.poisson(1, n_target),
    'referral_score': np.random.beta(1, 9, n_target),
    'treatment': np.random.binomial(1, 0.34, n_target),
    'date': pd.date_range('2024-01-01', periods=n_target, freq='D')
})

# 生成转化结果（真实ATE=3.2%）
true_ate = 0.032
target_data['converted'] = (
    0.08 +  # 基线转化率较低
    0.001 * target_data['browsing_pages'] +
    0.025 * target_data['referral_score'] +
    true_ate * target_data['treatment'] +
    np.random.normal(0, 0.05, n_target)
)
target_data['converted'] = (target_data['converted'] > 0).astype(int)

# 执行完整流程
system = CrossBorderCausalSystem(source_data, target_data)

# 阶段1-5顺序执行
print("=== 阶段1: 数据增强 ===")
system.phase1_data_augmentation()

print("\n=== 阶段2: 迁移学习 ===")
system.phase2_transfer_learning()

print("\n=== 阶段3: 贝叶斯建模 ===")
system.phase3_bayesian_modeling()

print("\n=== 阶段4: 合成控制验证 ===")
system.phase4_synthetic_control_validation()

# 预留池用于主动学习
candidate_pool = target_data[target_data['treatment'] == 0].copy()
print("\n=== 阶段5: 主动学习建议 ===")
system.phase5_active_learning_suggestion(candidate_pool, n_suggest=50)

# 生成报告
final_report = system.generate_report()

# 保存模型和报告
import joblib
joblib.dump(system, 'crossborder_causal_system.pkl')
with open('causal_report.txt', 'w') as f:
    f.write(final_report)

print("\n系统已保存至 crossborder_causal_system.pkl")