使用蒙特卡洛方法对几何分布进行序贯概率比检验 (SPRT)

Question

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

在 2 个复合假设的情况下，用计算机实施广义顺序测试，
使用蒙特卡洛方法估计其错误类型 I 和 II 的概率，以及条件预期样本量。

实现：几何分布，$\theta$ 是它唯一的参数。

分析错误概率和“假设接近度”和“大小”的预期样本大小的依赖关系。
将 $\alpha$ 和 $\beta$ 用于您选择的阈值计算。

我的解答思路和尝试过的方法，不写自己思路的，回答率下降 60%

我找到的方法代码如下，但是仅仅只是实现几何分布的蒙特卡洛方法，不知道该怎样实现使用蒙特卡洛方法对几何分布进行序贯概率比检验 (SPRT)。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(222)

# 把计算得到的函数写成一个函数
def distribution_z(z, p, max_k=200):
    import math
    j = int(math.floor(z))
    A = 0
    for m in range(1, j + 1):
        A += (1 - p) ** (m - 1)
    A *= p

    B = 0
    for k in range(j + 1, max_k + 1):
        a = (1 - p) ** (k - 1)
        a /= k
        B += a
    B *= z * p

    return A + B


def pdf_z(z, p, max_k=200):
    import math
    j = int(math.floor(z))
    B = 0
    for k in range(j + 1, max_k + 1):
        a = (1 - p) ** (k - 1)
        a /= k
        B += a
    return B * p


p = 0.1
# 选取数据点，点越多越精确
dataPoints = 10000

Unit = np.random.rand(dataPoints)
Geom = np.random.geometric(p, dataPoints)
distri_of_Monte = Geom * Unit

# 概率密度函数 PDF
plt.hist(distri_of_Monte, bins=40, range=(0, 40))
points_of_z = np.arange(0, 41, 0.01)
pdf_of_z = np.array([pdf_z(zi, p) for zi in points_of_z]) * dataPoints
plt.plot(points_of_z, pdf_of_z)
# print(pdf_of_z)
plt.show()

hist, bin_edges = np.histogram(distri_of_Monte, bins=40, range=(0, 40))

# 概率分布函数 CDF
hist_list = np.cumsum(hist) / dataPoints

plt.plot(bin_edges[1:], hist_list)

points_of_z = np.arange(1, 41, 0.1)
distri_of_z = [distribution_z(zi, p) for zi in points_of_z]

plt.plot(points_of_z, distri_of_z)

plt.show()

import sprt as sprt
import numpy as np

# Null value
h0 = 0.5
# Alternative value
h1 = 0.55
# Type I error rate = 0.05
alpha = 0.05
# Type II error rate = 0.2
beta = 0.2
# Values
values = np.random.binomial(1, 0.55, 100)
test = sprt.SPRTBinomial(h0 = h0, h1 = h1, alpha = alpha, beta = beta, values = values)

test.plot()

# Plot the data and boundary but without fill the color
test.plot(fill = False)

我想要达到的结果，如果你需要快速回答，请尝试 “付费悬赏”

使用蒙特卡洛方法，对几何分布进行序贯概率比检验 (SPRT)，同时满足题目的要求。使用Python编写代码。

Answer 1

可以先导入所需的库：

import numpy as np
import matplotlib.pyplot as plt

然后可以定义一个函数来实现序贯概率比检验 (SPRT) 算法。这个函数接受一组观察数据和两个阈值 $\alpha$ 和 $\beta$ 作为参数，并返回是否应该拒绝假设。

def sprt(data, alpha, beta):
    n = len(data)
    p_hat = data.mean()
    z = p_hat / (1 - p_hat)
    a = (1 - beta) / alpha
    b = beta / (1 - alpha)
    for i in range(n):
        z *= data[i] / (1 - data[i])
        if z < a:
            return "Reject H0"
        elif z > b:
            return "Reject H1"
    return "Fail to reject"

然后可以使用蒙特卡罗方法来模拟数据生成过程。假设已经实现了一个函数来计算几何分布的概率密度函数（pdf）和累积分布函数（cdf），这里就不再赘述。

接下来，可以定义一个函数来生成符合几何分布的观察数据。这个函数接受几何分布的唯一参数 $\theta$ 和要生成的数据点数量 $n$ 作为参数，并返回一组符合几何分布的观察数据。

def generate_data(theta, n):
    p = 1 - theta
    data = np.random.geometric(p, n)
    return data

然后可以使用这个函数来生成大量的符合几何分布的观察数据，并使用 SPRT 算法对这些数据进行检验。可以记录每次检验的结果，并统计错误类型 I 和 II 的次数。

# 设置参数
theta = 0.5
alpha = 0.05
beta = 0.1
n_trials = 10000

# 记录错误类型 I 和 II 的次数
error_I = 0
error_II = 0

# 运行多次模拟
for i in range(n_trials):
    # 生成符合几何分布的观察数据
    data = generate_data(theta, n)
    # 进行序贯概率比检验
    result = sprt(data, alpha, beta)
    if result == "Reject H0":
        error_I += 1
    elif result == "Reject H1":
        error_II += 1

# 计算错误概率
error_I_prob = error_I / n_trials
error_II_prob = error_II / n_trials

可以使用计算出的错误概率来估计条件预期样本大小。可以通过计算错误概率与假设接近度之间的依赖关系来估计条件预期样本大小。这里假设已经定义了一个函数 expected_sample_size 来计算条件预期样本大小。

# 计算条件预期样本大小
n_0 = expected_sample_size(error_I_prob, alpha)
n_1 = expected_sample_size(error_II_prob, beta)

望采纳。