如何用scipy获得指数分布的置信区间？

Question

在指数分布中使用 scipy 获取置信区间的正确方法是什么？

from scipy.stats import expon

conf_int = scipy.stats.expon.ttfs(0.95, loc=0, scale=1)

Answer 1

让我们澄清一下 scipy.stats 的一般用法

1. 如果您说

   from scipy.stats import expon

   ，则只有

   expon

   类将在您的命名空间中作为

   expon

   使用。

2. 如果您想使用

   scipy.stats.expon

   ，则必须使用

   import scipy.stats

   ，这样您就可以使用

   scipy.stats

   ，并且从那里您可以访问它的任何类。

3. expon

   是一个 rv_连续，它是实值连续分布的通用接口。所以你应该坚持这些方法，显然

   ttfs

   不可用。

4. 假设你想要一个居中的 95% 区间，那么你需要的是分位数

   2.5%

   和分位数

   100%-2.5%=97.5%

5. 您可以使用百分点函数 ppf 获得任何分位数，将百分比表示为

   100%

   的分数

将它们放在一起，您最终应该得到如下代码

from scipy.stats import expon
expon.ppf([0.025, 0.975], loc=0, scale=1)

array([0.02531781, 3.68887945])

您还可以构造一个代表参数化分布的

expon

实例

distribution = expon(loc=0, scale=1)
distribution.ppf([0.025, 0.975])

或者为系数区间选择不同的位置，例如

expon.ppf([0.0, 0.95], loc=0, scale=1)

array([0.        , 2.99573227])

Answer 2

指数分布对于速率参数 lambda 有一个封闭式置信区间

 公式 

。我认为考虑 CI 与您关心的参数相关始终很重要，或者您只是为了描述目的而查找分位数值。

使用Python，你可以做这样的事情：

import numpy as np
from scipy.stats import chi2
from typing import Tuple

def ci_exp(sample: np.ndarray, alpha=0.05) -> Tuple[float, float]:
    """
    Finds the confidence interval for the rate parameter of an exponential  distribution.

    Parameters
    ----------
    samples : np.ndarray
        The sample of realized exponential random variables.
    alpha : float, optional
        The level of significance, by default 0.05.

    Returns
    -------
    Tuple[float, float]
        The lower and upper bounds of the confidence interval.
    """
    n = len(sample)
    # MLE of lambda
    sample_mean = np.mean(sample)

    # Compute the sum of samples
    sample_sum = np.sum(sample)

    # The sample sum is gamma distributed
    gamma_shape = n
    gamma_scale = 1 / sample_mean

    chi2_df = 2 * n

    # Confidence interval for lambda
    lower_bound = (2 * sample_sum) / chi2.ppf(q=(1 - (alpha / 2)), df=chi2_df)
    upper_bound = (2 * sample_sum) / chi2.ppf(q=alpha / 2, df=chi2_df)
    ci = (lower_bound, upper_bound)

    return ci

if __name__ == "__main__":
    sample = np.random.exponential(scale=1/0.5, size=100)
    ci = ci_exp(sample)
    print(ci)