基于1/2-exp(1)误差的加速度系数对称置信区间优化问询

工厂针对某机型开展加速度系数测试：选取n台该机型，测量其38秒内的行驶距离。已知路径测量误差独立同分布，服从1/2-exp(1)分布，需要构建加速度系数的对称置信区间。输入参数为置信水平（0-1之间的数值）与路径测量值的一维numpy数组，返回值为置信区间左右边界组成的元组或列表。

现有Python实现代码得到的区间长度过长，需进行优化：

import pandas as pd
import numpy as np
from scipy.stats import norm

def solution(p: float, x: np.array) -> tuple:

    x_bar = np.mean(x)
    s = np.std(x, ddof=1)
    z = norm.ppf(1 - (1 - p)/2)
    lower = x_bar - z * s / np.sqrt(len(x))
    upper = x_bar + z * s / np.sqrt(len(x))
    t = 38
    v = x.mean() / t
    a = 2 * v / t
    lower_a = 2 * lower / t
    upper_a = 2 * upper / t
    return lower_a, upper_a

问题根源与优化思路 #

• 误用正态分布近似：现有代码直接使用正态分布Z分位数构建区间，但测量误差服从的是1/2指数分布（均值为1/2、方差为1/4），并非正态分布。当样本量较小时，正态近似会过度放大区间宽度。
• 利用误差分布的精确性质：
  设真实行驶距离为( S_i )，测量值( X_i = S_i + \epsilon_i )，其中( \epsilon_i \sim \text{Exp}(1)/2 )，即( 2\epsilon_i \sim \text{Exp}(1) )。
  结合加速度公式( S_i = \frac{1}{2} a t^2 )，可推导出( 2X_i - a t^2 \sim \text{Exp}(1) )。n个该变量之和( T = 2n\bar{X} - n a t^2 )服从Gamma分布( \text{Gamma}(n, 1) )。
• 构建精确置信区间：
  针对置信水平( p )，找到Gamma(n,1)分布的( (1-p)/2 )和( 1-(1-p)/2 )分位数( q_{\text{low}} )、( q_{\text{high}} )，通过不等式( P(q_{\text{low}} \leq T \leq q_{\text{high}}) = p )反推得到( a )的置信区间。

优化后代码 #

import numpy as np
from scipy.stats import gamma

def solution(p: float, x: np.array) -> tuple:
    n = len(x)
    x_bar = np.mean(x)
    t = 38
    # 获取Gamma(n,1)分布的双侧分位数
    q_low = gamma.ppf((1 - p)/2, a=n, scale=1)
    q_high = gamma.ppf(1 - (1 - p)/2, a=n, scale=1)
    # 加速度系数的点估计
    a_hat = 2 * x_bar / (t ** 2)
    # 计算置信区间边界
    lower_a = a_hat - q_high / (n * t ** 2)
    upper_a = a_hat - q_low / (n * t ** 2)
    return lower_a, upper_a

内容的提问来源于stack exchange，提问作者Joker221