两参数指数分布拟合单参数指数分布的LRT中Type I Error Rate高于显著性水平的原因排查

看起来你遇到的问题主要源于LRT经典渐近卡方假设的失效，再加上代码中的一个小错误（虽未影响最终拟合结果，但逻辑上需要修正）。让我一步步拆解原因和解决方案：

1. 核心原因：LRT的正则条件不满足，渐近卡方分布假设完全失效 #

经典LRT统计量渐近服从卡方分布的前提是H0下的真实参数位于全参数空间的内部，且模型满足一系列正则条件（比如Fisher信息矩阵非奇异）。但在你的问题中：

• 两参数指数分布的对数似然对位置参数θ的一阶导数是常数（$n\lambda$），这直接导致Fisher信息矩阵中θ对应的方差项为0，矩阵奇异——这违反了正则条件的核心要求。
• 当H0成立（θ=0）时，两参数指数分布的θ的MLE是样本最小值（而非真实值0），这进一步偏离了LRT渐近理论的假设场景：渐近卡方分布要求H0下的MLE收敛到真实参数，而这里的MLE却收敛到样本最小值（一个随样本变化的统计量）。

在这种情况下，LRT统计量的渐近分布绝对不是卡方(1)，而是完全不同的分布。你强行用卡方(1)计算p值，自然会导致Type I error率严重偏离预期。

2. H0下统计量的真实分布与卡方(1)差异极大 #

我们可以推导H0下你的LRT统计量的具体形式：
对于来自单参数指数分布（θ=0）的样本，设样本均值为$\mu$，最小值为$m$，统计量：
$$\Lambda = 2n\log\left(\frac{\mu}{\mu - m}\right)$$
当样本量n=100时，这个统计量的分布和卡方(1)差异极大：卡方(1)的95%分位数是2.706，但你的统计量有远高于5%的概率超过这个值——这就是为什么你观察到Type I error率达到0.15（远高于设定的0.05）：很多在真实分布中不显著的统计量，被卡方(1)的临界值误判为显著。

3. 代码中的参数边界错误（逻辑问题，未影响拟合结果） #

你在拟合全模型时，将theta的边界设置为(np.min(data), np.max(data))，这要求θ≥样本最小值，但两参数指数分布的参数空间是θ≤样本最小值（所有数据点必须≥θ）。
不过这个错误没有改变拟合结果：当θ>样本最小值时，你的log似然函数会返回-无穷，优化器只能被迫选择θ=样本最小值（这正是两参数指数θ的MLE）。但为了代码逻辑的正确性，建议修正为：

bounds=[(-np.inf, np.min(data)), (0, None)]

如何修正Type I error率偏高的问题 #

• 使用经验分布计算p值：针对给定的样本量n，模拟大量H0下的数据集，计算Lambda统计量并得到其经验分布。实际检验时，用这个经验分布来计算p值（比如，统计实际Lambda在经验分布中的分位数）。示例代码如下：

  # 预先模拟H0下Lambda统计量的经验分布
  def simulate_null_distribution(sample_size, lambda_true, num_sim=10000):
      lambdas = []
      for _ in range(num_sim):
          data = np.random.exponential(1/lambda_true, sample_size)
          m = np.min(data)
          mu = np.mean(data)
          # 计算全模型和简化模型的对数似然
          ll_full = -sample_size * np.log(mu - m) - sample_size
          ll_reduced = -sample_size * np.log(mu) - sample_size
          lam = 2 * (ll_full - ll_reduced)
          lambdas.append(lam)
      return np.array(lambdas)
  
  # 使用经验分布计算p值
  null_dist = simulate_null_distribution(sample_size, lambda_true)
  def likelihood_ratio_test_empirical(data, null_dist):
      m = np.min(data)
      mu = np.mean(data)
      ll_full = -len(data) * np.log(mu - m) - len(data)
      ll_reduced = -len(data) * np.log(mu) - len(data)
      lam = 2 * (ll_full - ll_reduced)
      # 计算p值：经验分布中大于等于当前Lambda的比例
      p_value = np.mean(null_dist >= lam)
      return p_value
  

• 考虑替代检验方法：比如Score检验或Wald检验，但同样需要注意它们的渐近假设是否成立；或者直接使用针对该特定问题的精确检验方法。

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