很高兴看到你已经把两阶段Bootstrap的实现落地了，咱们逐个拆解你的问题：

1. 两阶段Bootstrap方法是否合理？ #

这个思路完全站得住脚！你用的是非参数+参数混合的Bootstrap，两个阶段的分工非常清晰：

• 第一阶段的非参数重采样（按索引有放回抽样）：用来捕捉原始样本的统计变异性——也就是如果你的数据是从某个总体中抽出来的子集，这一步能模拟你再抽一次样本会得到什么结果。
• 第二阶段的参数抽样（从$\mathcal{N}(y_i, e_i^2)$生成合成数据）：用来还原每个数据点本身的测量不确定性——毕竟$y_i$只是真实值的带误差观测，不是真值本身。

这种方法完美兼顾了两种误差来源，尤其适合你既想保留原始样本的分布特征（非参数部分），又要严谨考虑测量误差的场景，逻辑上是完全通顺的。

2. 权重选择：$w_i=1/e_i$还是$w_i=1/e_i^2$？ #

这里要纠正一下：正确的权重应该是$w_i=1/e_i^2$（逆方差权重），而不是逆标准差。

原因来自高斯假设下的最优估计（Gauss-Markov定理）：当每个观测的测量误差方差为$e_i^2$时，逆方差权重能让加权均值成为总体均值的最小方差无偏估计（MVUE）。举个例子：如果A点的误差是B点的2倍，那A点的方差是B点的4倍，它对均值的贡献应该是B点的1/4，而不是1/2——精度越低的点，权重应该按方差的倒数衰减，这样才能最小化最终均值的方差。

你只需要把代码里的权重计算改成：

w_draw = 1.0 / (e_draw ** 2)

对应的加权均值计算保持sum(y_synth * w_draw) / sum(w_draw)就可以了。

3. 有没有更标准或高效的替代方法？ #

根据你的需求（同时考虑样本抽样变异性+测量误差），你的Bootstrap方法已经是很合适的选择，但可以补充几种场景化的替代方案：

（1）仅考虑测量误差：加权最小二乘（WLS）直接估计

如果你的原始样本就是全部研究对象（没有抽样误差，只有测量误差），那可以跳过Bootstrap，直接用WLS估计：

• 均值估计：$\hat{\mu} = \frac{\sum (y_i / e_i^2)}{\sum (1 / e_i^2)}$
• 标准误差：$\text{SE}(\hat{\mu}) = \frac{1}{\sqrt{\sum (1 / e_i^2)}}$

这个方法计算极快，但只适用于“无样本抽样误差”的场景，无法覆盖你需要考虑统计变异性的需求。

（2）参数化替代：随机效应模型

如果你的数据符合“每个数据点的真实值围绕总体均值服从高斯分布”的假设（即真实值$\mu_i \sim \mathcal{N}(\mu, \tau^2)$，观测值$y_i \sim \mathcal{N}(\mu_i, e_i^2)$），可以用元分析中的随机效应模型，通过极大似然估计直接得到$\mu$的估计值和置信区间。

但这个方法依赖于高斯假设的合理性，而你的Bootstrap是非参数的，不需要任何分布假设，灵活性更强。

（3）简化版Bootstrap（当无抽样误差时）

如果确认不需要考虑样本的统计变异性，可以去掉第一阶段的非参数重采样，直接对每个$y_i$重复从$\mathcal{N}(y_i, e_i^2)$抽样，再计算加权均值——这会大幅提升计算效率，因为少了一次重采样步骤。

代码修改建议（关键部分） #

把bootstrap_iteration里的权重计算改为逆方差：

def bootstrap_iteration(_):
    # Nonparametric resample
    indices = rng.integers(0, N, size=N)
    y_draw  = y_vals[indices]
    e_draw  = e_vals[indices]

    # Parametric resample from each point
    y_synth = rng.normal(loc=y_draw, scale=e_draw)

    # 修正为逆方差权重
    w_draw  = 1.0 / (e_draw ** 2)
    w_sum   = np.sum(w_draw)
    if w_sum == 0:
        return np.nan
    return np.sum(y_synth * w_draw) / w_sum

备注：内容来源于stack exchange，提问作者Jokerp