实施CUSUM求解规则，以检测观察流中的不连续性

Question

问题遇到的现象和发生背景

实施CUSUM求解规则，以检测观察流中的不连续性（观察值的概率分布中的一次性跳跃变化）。观察流的数据应该用 "没有不连续 "和 "在时间$/nu$有不连续 "两种方式建模。
与提问问题1和问题2类似，使用蒙特卡洛方法为已实现的求解器估计误报概率和故障后到做出决定的平均延迟（以观测次数计）。考虑在没有篡改的情况下对所构建的解算器的约束（避免无限等待）。
独立设置问题的参数：阈值、衰减时间等。
调查阈值和两个备选参数值的接近程度对性能的影响（错误警报的概率，故障后决定前的平均失踪观测值数量）。
问题1
在 2 个简单假设的情况下，用计算机实施序贯概率比检验 (SPRT)。
使用蒙特卡洛方法估计其错误类型 I 和 II 的概率，以及条件预期样本大小（运行测试 10 000 次，$H_0$ 为真，
计算有利于 $H_1$ 的 FALSE 决策的数量，除以 10 000 -> 你得到 $\alpha$ 的估计值；启动一个计数器，
将你停止测试的观察次数添加到一个计数器，每次运行时 $H_0$ 为真，然后除以 10 000 ->
你将得到 $ 的估计值t_0$；$\beta$ 和 $t_1$ 的估计是相似的）。

实现：几何分布，$\theta$ 是它唯一的参数。
分析错误概率和假设接近度的预期样本大小的依赖关系。将 $\alpha$ 和 $\beta$ 用于您选择的阈值计算

尝试过的解决方法

与上两个提问问题类似

我想要达到的结果

用Python编写写一下注释，使用蒙特卡洛方法，对几何分布进行检验，同时满足题目的要求。

Answer 1

基于Monster 组和GPT的调写：
用Python实现的序贯概率比检验的示例，针对几何分布的情况。代码中使用蒙特卡罗方法估计错误类型 I 和 II 的概率，以及条件预期样本大小。

import numpy as np

def sequential_probability_ratio_test(sample, threshold, alpha, beta):
    # 初始化统计量
    t = 0
    n = len(sample)
    log_l = 0
    log_f = 0
    decision = 0

    # 计算参数
    p0 = 1 / (1 + threshold)
    p1 = threshold / (1 + threshold)

    # 进行SPRT测试
    while decision == 0 and t < n:
        x = sample[t]
        log_l += np.log(p1 * (1 - p0)**x)
        log_f += np.log(p0 * (1 - p1)**x)
        t += 1

        # 判断是否需要停止测试
        if log_l - log_f >= np.log(beta / alpha):
            decision = 1
        elif log_f - log_l >= np.log(beta / alpha):
            decision = -1

    # 返回统计信息
    return decision, t

# 进行蒙特卡罗模拟
n_simulations = 10000
alpha = 0.05
beta = 0.1
theta = 0.2
threshold = 2

n_false = 0
n_miss = 0
false_count = 0
miss_count = 0

for i in range(n_simulations):
    sample_size = sequential_probability_ratio_test(np.random.geometric(theta, 100), threshold, alpha, beta)[1]
    if sample_size < 100:
        n_false += 1
        false_count += 100 - sample_size
    else:
        n_miss += 1
        miss_count += sample_size - 100

# 计算误差率和平均停顿时间
false_rate = n_false / n_simulations
miss_rate = n_miss / n_simulations
avg_false_count = false_count / n_false if n_false > 0 else 0
avg_miss_count = miss_count / n_miss if n_miss > 0 else 0

# 打印结果
print("错误类型 I 的概率：", false_rate)
print("错误类型 II 的概率：", miss_rate)
print("条件预期样本大小（平均失踪观测值数量）：", avg_miss_count)
print("错误警报的数量：", avg_false_count)

# 计算阈值为2时的 alpha 和 beta
p0 = 1 / (1 + threshold)
p1 = threshold / (1 + threshold)
a = (np.log(beta / alpha) + np.log((1 - p1) / p1)) / np.log((1 - p0) / p0)
b = (np.log(alpha / beta) + np.log((1 - p0) / p0)) / np.log((1 - p1) / p1)
print("阈值为2时的 alpha 和 beta：", alpha**a, beta**b)

首先定义了一个名为 sequential_probability_ratio_test 的函数，该函数接受样本、阈值、α和β等参数，然后进行SPRT测试并返回决策和观测样本数等统计信息。

然后进行蒙特卡罗模拟，使用 np.random.geometric 生成一个参数为 theta 的几何分布样本，并调用 sequential_probability_ratio_test 函数进行SPRT测试。如果测试结果表明样本数小于100，则认为是错误类型 I，统计错误数和错误警报的数量；否则认为是错误类型 II，统计错失数和条件预期样本大小等信息。

最后，计算错误类型 I 和 II 的概率，以及条件预期样本大小和错误警报的数量等统计信息，并打印输出。

此外，该代码还计算了阈值为2时的α和β，并将其打印输出。