当前位置：首页 > article >正文

使用 scipy 计算置信区间

article 2024/11/14 7:43:15

问题

我有一个形状为的数组(n, timesteps)，其中n是试验次数，timesteps是每次试验的长度。该数组的每个值表示一个随机测量值。
我想实现一个通用函数，该函数计算给定统计数据（平均值、中位数……）的置信区间，假设 1）基础分布为正态分布，2）为学生分布。

类似于

def normal_ci(
    data: np.array,
    axis: int = 0,
    statistic: Callable = np.mean,
    confidence: float = 0.95
):

以及类似的功能student_ci()。

我的问题是，默认情况下，scipy 函数会计算平均值的区间，对吗？

解决方案

计算答案是使用bootstrap。

from typing import Callable
import numpy as np
from scipy import stats

rng = np.random.default_rng(84354894298246)
data = rng.standard_normal(size=(1000, 100))

def normal_ci(
    data: np.array,
    axis: int = 0,
    statistic: Callable = np.mean,
    confidence: float = 0.95
):
    res = stats.bootstrap((data,), statistic, axis=axis,
                          confidence_level=confidence)
    return tuple(res.confidence_interval)

low, high = normal_ci(data, axis=-1)
np.mean((low < 0) & (0 < high))  # 0.953

由于您知道数据来自哪些家族，因此您可以研究参数引导法。

def normal_ci(
    data: np.array,
    axis: int = 0,
    statistic: Callable = np.mean,
    confidence: float = 0.95
):
    # fit a normal distribution to the data
    mean = np.mean(data, axis=axis)
    std = np.std(data, ddof=1, axis=axis)

    # resample data from the fitted normal distribution
    n_resamples = 999
    m = data.shape[axis]  # assuming axis is an integer
    resample = rng.normal(loc=mean, scale=std, size=(m, n_resamples) + mean.shape)

    # compute the statistic for each of the resamples to estimate
    # the distribution
    statistic_distribution = statistic(resample, axis=0)

    # Generate the confidence interval
    # percentile bootstrap (very crude)
    alpha = 1 - confidence
    low, high = np.quantile(statistic_distribution, [alpha/2, 1-alpha/2], axis=0)
    return low, high

low, high = normal_ci(data, axis=-1)
np.mean((low < 0) & (0 < high))  # 0.954