2025-11-05-基于随机加权推断模型（IM）的参数估计算法

基于随机加权推断模型（IM）的参数估计方法#

不可观测的泊松过程 $B \sim \text{Poisson}(b)$ 和 $S \sim \text{Poisson}(q)$ ，观测值 $Y = B + S \sim \text{Poisson}(b + q)$ 和辅助变量 $W \sim \text{Poisson}(m \cdot b)$ ，目标是估计参数 $b$ 、 $q$ m?

现提出基于随机加权推断模型（Inferential Model, IM）的算法通过引入随机权重处理离散分布，来构造置信区间

方法概述#

随机加权 IM 方法的核心是通过引入均匀随机变量和权重，将离散泊松分布连续化，从而逆解参数该方法假设 $m$ 已知算法分为四个步骤：建立关联模型、引入随机权重连续化、逆解参数、模拟验证覆盖率

算法步骤#

1: 建立关联模型#

对观测值 $Y$ ，其分布函数为 $F_{θ}(y)$ ，其中 $θ = b + q$ 引入辅助变量 $u \sim \text{Uniform}(0,1)$ ，建立关联：

$F_{θ}(Y-1) \leq u < F_{θ}(Y)$

对观测值 $W$ ，其分布函数为 $F_{mb}(w)$ ，引入辅助变量 $v \sim \text{Uniform}(0,1)$ ，建立关联：

$F_{mb}(W-1) \leq v \leq F_{mb}(W)$

2: 引入随机权重进行连续化#

引入随机权重 $w_1, w_2 \sim \text{Uniform}(0,1)$ ， $u, v \sim \text{Uniform}(0,1)$ 构造方程：

$G(θ): w_1 F_{θ}(Y-1) + (1-w_1) F_{θ}(Y) = u$

$H(mb): w_2 F_{mb}(W-1) + (1-w_2) F_{mb}(W) = v$

$G$ 和 $H$ 是关于参数的函数，且由于泊松分布函的单调性，逆函数 $G^{-1}$ 和 $H^{-1}$ 存在

3: 逆解参数，使用 brentq 函数进行二分法求根#

从上述方程解出 $θ$ 和 $mb$ ：

$θ = G^{-1}(u), \quad mb = H^{-1}(v)$

推导参数：

$b = \frac{H^{-1}(v)}{m}, \quad q = θ - b = G^{-1}(u) - \frac{H^{-1}(v)}{m}$

得到 $b$ 和 $q$ 的表达式

4: 验证覆盖率#

模拟目标：检查构造的置信区间覆盖真实 $q$ 的概率是否接近 95%
步骤：

生成观测数据：从分布生成一对观测值 $(Y,W)$

对于固定 $(Y,W)$ ，生成大量样本（如 $N=10000$ 个）的 $u, v, w_1, w_2 \sim \text{Uniform}(0,1)$

计算 $q$ 的候选值：对于每个样本，计算 $q_{\text{candidate}} = G^{-1}(u) - \frac{H^{-1}(v)}{m}$

构造置信区间：从 10000 个 q 值中取分位数，得到置信区间 $[q\_{0.025}, q\_{0.975}]$

检查覆盖率：重复步骤 2-5 多次（ $M=10000$ 次），每次生成新的 $(Y,W)$ ，计算区间覆盖真实 $q$ 的比例覆盖率 $= \frac{\text{覆盖次数}}{M}$

算法实现#

1
import numpy as np
2
import pandas as pd
3
from scipy.stats import poisson, gamma
4
from scipy.optimize import brentq
5
import matplotlib.pyplot as plt
6
import seaborn as sns
7
from tqdm import tqdm
8
import warnings
9
warnings.filterwarnings('ignore')
10

11

12
# 设置全局字体为 SimHei (黑体) 或其他中文字体
13
plt.rcParams['font.sans-serif'] = ['SimHei']
14
plt.rcParams['axes.unicode_minus'] = False
15

16
class PoissonSignalInference:
17
    """
18
    基于随机加权IM方法的泊松信号推断类
19
    用于估计背景率b、信号率q和尺度参数m
20
    """
21

22
    def __init__(self, b_true=None, q_true=None, m_true=None):
23
        """
24
        初始化真实参数值
25
        """
26
        self.b_true = b_true
27
        self.q_true = q_true
28
        self.m_true = m_true
29
        self.STA_true = b_true + q_true if b_true is not None and q_true is not None else None
30

31
    def generate_data(self, n_samples=1):
32
        """
33

34
        生成观测数据Y和W
35

36
        当前有不可观测的泊松过程 B ∼ Poisson(b) 和 S ∼ Poisson(q)
37
        可观测对象有 Y = B + S ∼ Poisson(b + q) 和辅助变量 W ∼ Poisson(m · b)
38
        目标是估计参数 b、q 和可能的 m
39

40
        """
41
        if self.b_true is None or self.q_true is None or self.m_true is None:
42
            raise ValueError("必须先设置真实参数值")
43

44
        if self.STA_true is None:
45
            self.STA_true = self.b_true + self.q_true
46

47
        Y = np.random.poisson(self.STA_true, n_samples)
48
        W = np.random.poisson(self.m_true * self.b_true, n_samples)
49

50
        return Y, W
51

52
    def poisson_cdf(self, k, theta):
53
        """
54
        计算泊松分布的累积分布函数
55
        """
56
        return poisson.cdf(k, theta) #泊松分布的CDF公式，参数为 k 和 theta
57

58
    def solve_G_inverse(self, u, Y, w1, bracket_low=0, bracket_high=100):
59
        """
60
        求解方程 G(STA) = u 的逆函数，得到STA
61
        使用二分法求解 G(STA): w_1 F_{STA}(Y-1) + (1-w_1) F_{STA}(Y) = u
62

63
        计算泊松混合分布的累积分布函数值与目标值的差值
64

65
        该函数计算一个泊松混合分布的CDF值，该混合分布由两个泊松分布组成，权重分别为w1和(1-w1)，
66
        然后减去目标值u，通常用于求解分位数或进行概率计算
67
        参数:
68
            Y (int or array-like): 观测值，泊松分布的计数变量
69
            STA (float or array-like): 泊松分布的参数λ（均值）
70
            w1 (float): 第一个泊松分布的权重，取值范围[0,1]
71
            u (float): 目标概率值，用于比较或求解
72
        返回:
73
        """
74
        def equation(STA): #返回值:(float or array-like): 泊松混合分布的CDF值与目标值u的差值
75
            if Y == 0:
76
                # 当Y=0时，F_STA(-1)=0，方程简化为F_STA(0)=u
77
                return self.poisson_cdf(0, STA) - u
78
            else:
79
                return w1 * self.poisson_cdf(Y-1, STA) + (1-w1) * self.poisson_cdf(Y, STA) - u
80
                # 求加权泊松混合分布
81
                # 第一项：权重w1乘以P(X <= Y-1)的概率
82
                # 第二项：权重(1-w1)乘以P(X <= Y)的概率
83
                # 最后减去目标值u
84

85
        try:
86
            """
87
            使用二分法求解 G(STA) = u 方程中的 STA 值
88
            """
89
            # 设置上下界来初始化搜索区间
90
            low = max(0.00001, bracket_low) #bracket_low下界（）
91
            high = max(10 * Y, bracket_high) if Y > 0 else bracket_high #bracket_high上界，根据观测值 Y 调整上界
92

93
            # 检查边界值
94
            f_low = equation(low)
95
            f_high = equation(high)
96

97

98
            if f_low * f_high > 0:
99
                # 如果 f_low * f_high > 0（同号），则自动扩大搜索区间
100
                attempts = 0
101
                while f_low * f_high > 0 and attempts < 10:
102
                    high *= 2 #最多尝试10次，每次将上界翻倍：high *= 2
103
                    f_high = equation(high)
104
                    attempts += 1
105

106
            STA_solution = brentq(equation, low, high, xtol=1e-6)
107
            """
108
            参数为要求解的方程（加权泊松混合分布），上下区间，精度
109
            使用 brentq 函数进行二分法求根
110
            求解 w1 * self.poisson_cdf(Y-1, STA) + (1-w1) * self.poisson_cdf(Y, STA) - u = 0
111

112
            即求解方程 G(STA) = u
113
            得到STA
114

115
            """
116

117
            return STA_solution #返回求解得到的 STA 值
118
        except (ValueError, RuntimeError):
119
            # 如果求解失败，返回保守估计
120
            print("ERR solve_G_inverse 求解失败，返回保守估计")
121
            return max(0.001, Y)
122

123
    def solve_H_inverse(self, v, W, w2, m, bracket_low=0, bracket_high=100):
124
        """
125
        求解方程 H(mb) = v 的逆函数，得到mb
126

127
        H(mb): w_2 F_{mb}(W-1) + (1-w_2) F_{mb}(W) = v
128

129
        """
130
        def equation(mb):
131
            if W == 0:
132
                return self.poisson_cdf(0, mb) - v
133
            else:
134
                return w2 * self.poisson_cdf(W-1, mb) + (1-w2) * self.poisson_cdf(W, mb) - v
135

136
        try:
137
            low = max(0.001, bracket_low)
138
            high = max(10 * W, bracket_high) if W > 0 else bracket_high
139

140
            f_low = equation(low)
141
            f_high = equation(high)
142

143
            if f_low * f_high > 0:
144
                attempts = 0
145
                while f_low * f_high > 0 and attempts < 10:
146
                    high *= 2
147
                    f_high = equation(high)
148
                    attempts += 1
149

150
            mb_solution = brentq(equation, low, high, xtol=1e-6)
151
            return mb_solution
152
        except (ValueError, RuntimeError):
153
            return max(0.001, W)
154

155
    def estimate_parameters_IM(self, Y_obs, W_obs, m_known=True, n_inner_samples=10000):
156
        """
157
        使用随机加权IM方法估计三个参数q、b和STA
158
        Y_obs：观测值Y
159
        W_obs：观测值W
160
        m_known：一个布尔值，表示参数m是否已知
161
        n_inner_samples：内部采样次数，默认为10000
162

163
        """
164
        if not m_known:
165
            raise NotImplementedError("m未知")
166

167
        # 创建三个列表存储候选参数值
168
        q_candidates = []
169
        b_candidates = []
170
        STA_candidates = []
171

172
        for _ in range(n_inner_samples):
173
            # 生成4个均匀分布的随机权重（u, v, w1, w2）
174
            u = np.random.uniform(0, 1)
175
            v = np.random.uniform(0, 1)
176
            w1 = np.random.uniform(0, 1)
177
            w2 = np.random.uniform(0, 1)
178

179
            # 使用solve_G_inverse和solve_H_inverse方法计算STA和mb的估计值
180
            STA_est = self.solve_G_inverse(u, Y_obs, w1)
181
            mb_est = self.solve_H_inverse(v, W_obs, w2, self.m_true)
182
            # G(STA): w_1 F_{STA}(Y-1) + (1-w_1) F_{STA}(Y) = u ==> 求解STA
183
            # H(mb): w_2 F_{mb}(W-1) + (1-w_2) F_{mb}(W) = v ==> 求解mb
184
            # 可以得到
185
            # STA = G⁻¹(u), mb = H⁻¹(v)
186
            # 计算b和q
187
            # b = H⁻¹(v)/m
188
            # q = STA - b = G⁻¹(u) - H⁻¹(v)/m
189
            if self.m_true is None:
190
                raise ValueError("参数 m_true 未设置，无法计算 b_est")
191
            b_est = mb_est / self.m_true
192
            q_est = STA_est - b_est
193

194
            # 确保非负
195
            q_est = max(0, q_est)
196
            b_est = max(0, b_est)
197

198
            #每次计算出新的 q_est 值后，会将其添加到 q_candidates 列表中
199
            q_candidates.append(q_est)
200
            b_candidates.append(b_est)
201
            STA_candidates.append(STA_est)
202

203
        return {
204
            'q_candidates': np.array(q_candidates),
205
            'b_candidates': np.array(b_candidates),
206
            'STA_candidates': np.array(STA_candidates)
207
        }
208

209
    def calculate_confidence_intervals(self, candidates_dict, alpha=0.05):
210
        """
211
        计算参数的置信区间
212
        """
213
        intervals = {}
214

215
        for param, values in candidates_dict.items():
216
            # 移除'_candidates'后缀，使键名简化
217
            simple_param = param.replace('_candidates', '')
218
            lower = np.quantile(values, alpha/2)
219
            upper = np.quantile(values, 1 - alpha/2)
220
            intervals[simple_param] = (lower, upper)
221

222
        return intervals
223

224
    def coverage_simulation(self, n_outer_trials=1000, n_inner_samples=1000, alpha=0.05):
225
        """
226
        进行覆盖率模拟验证
227
        """
228
        coverage_results = {
229
            'q': 0.0,
230
            'b': 0.0,
231
            'STA': 0.0
232
        }
233

234
        interval_lengths = {
235
            'q': [],
236
            'b': [],
237
            'STA': []
238
        }
239

240
        print("正在进行覆盖率模拟...")
241
        for i in tqdm(range(n_outer_trials)):
242
            # 生成新的观测数据
243
            Y_obs, W_obs = self.generate_data(1)
244
            Y_obs, W_obs = Y_obs[0], W_obs[0]
245

246
            # 估计参数
247
            candidates_dict = self.estimate_parameters_IM(Y_obs, W_obs,
248
                                                         n_inner_samples=n_inner_samples)
249

250
            # 计算置信区间
251
            intervals = self.calculate_confidence_intervals(candidates_dict, alpha)
252

253
            # 检查覆盖率
254
            if intervals['q'][0] <= self.q_true <= intervals['q'][1]:
255
                coverage_results['q'] += 1
256
            if intervals['b'][0] <= self.b_true <= intervals['b'][1]:
257
                coverage_results['b'] += 1
258
            if intervals['STA'][0] <= self.STA_true <= intervals['STA'][1]:
259
                coverage_results['STA'] += 1
260

261
            # 记录区间长度
262
            interval_lengths['q'].append(intervals['q'][1] - intervals['q'][0])
263
            interval_lengths['b'].append(intervals['b'][1] - intervals['b'][0])
264
            interval_lengths['STA'].append(intervals['STA'][1] - intervals['STA'][0])
265

266
        # 计算覆盖率
267
        for param in coverage_results:
268
            coverage_results[param] = coverage_results[param] / n_outer_trials
269

270
        return coverage_results, interval_lengths
271

272
    def plot_results(self, candidates_dict, intervals, Y_obs, W_obs):
273
        """
274
        可视化结果
275
        """
276
        fig, axes = plt.subplots(2, 2, figsize=(15, 10))
277

278
        # q的分布
279
        axes[0, 0].hist(candidates_dict['q_candidates'], bins=50, alpha=0.7, density=True)
280
        axes[0, 0].axvline(self.q_true, color='red', linestyle='--', linewidth=2, label=f'真实 q = {self.q_true}')
281
        axes[0, 0].axvline(intervals['q'][0], color='green', linestyle='--', linewidth=1, label=f'{100*(1-0.05)}% CI')
282
        axes[0, 0].axvline(intervals['q'][1], color='green', linestyle='--', linewidth=1)
283
        axes[0, 0].set_xlabel('q')
284
        axes[0, 0].set_ylabel('密度')
285
        axes[0, 0].set_title('q的后验分布')
286
        axes[0, 0].legend()
287

288
        # b的分布
289
        axes[0, 1].hist(candidates_dict['b_candidates'], bins=50, alpha=0.7, density=True)
290
        axes[0, 1].axvline(self.b_true, color='red', linestyle='--', linewidth=2, label=f'真实 b = {self.b_true}')
291
        axes[0, 1].axvline(intervals['b'][0], color='green', linestyle='--', linewidth=1, label=f'{100*(1-0.05)}% CI')
292
        axes[0, 1].axvline(intervals['b'][1], color='green', linestyle='--', linewidth=1)
293
        axes[0, 1].set_xlabel('b')
294
        axes[0, 1].set_ylabel('密度')
295
        axes[0, 1].set_title('b的后验分布')
296
        axes[0, 1].legend()
297

298
        # STA的分布
299
        axes[1, 0].hist(candidates_dict['STA_candidates'], bins=50, alpha=0.7, density=True)
300
        axes[1, 0].axvline(self.STA_true, color='red', linestyle='--', linewidth=2, label=f'真实 STA = {self.STA_true}')
301
        axes[1, 0].axvline(intervals['STA'][0], color='green', linestyle='--', linewidth=1, label=f'{100*(1-0.05)}% CI')
302
        axes[1, 0].axvline(intervals['STA'][1], color='green', linestyle='--', linewidth=1)
303
        axes[1, 0].set_xlabel('STA')
304
        axes[1, 0].set_ylabel('密度')
305
        axes[1, 0].set_title('STA的后验分布')
306
        axes[1, 0].legend()
307

308
        # 观测数据信息
309
        axes[1, 1].axis('off')
310
        text_str = f'观测数据:\nY = {Y_obs}\nW = {W_obs}\n\n真实参数:\nb = {self.b_true}\nq = {self.q_true}\nm = {self.m_true}\nSTA = {self.STA_true}'
311
        axes[1, 1].text(0.1, 0.9, text_str, transform=axes[1, 1].transAxes, fontsize=12,
312
                    verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
313

314
        plt.tight_layout()
315
        plt.show()
316

317
def main():
318
    """
319
    主函数：演示完整的推断流程
320
    """
321
    # 设置真实参数
322
    b_true = 3.0  # b
323
    q_true = 2.0  # q
324
    m_true = 0.8  # 尺度参数（已知）
325
    STA_true = b_true + q_true
326

327
    print("=== 泊松推断演示 ===")
328
    print(f"真实参数: b={b_true}, q={q_true}, m={m_true}, STA={STA_true}")
329

330
    # 初始化推断器
331
    inference = PoissonSignalInference(b_true, q_true, m_true)
332

333
    # 生成单次观测数据
334
    Y_obs, W_obs = inference.generate_data(1)
335
    Y_obs, W_obs = Y_obs[0], W_obs[0]
336
    print(f"\n生成的观测数据: Y={Y_obs}, W={W_obs}")
337

338
    # 使用IM方法进行参数估计
339
    print("\n正在进行参数估计...")
340
    candidates_dict = inference.estimate_parameters_IM(Y_obs, W_obs, n_inner_samples=10000)
341

342
    # 计算置信区间
343
    alpha = 0.05  # 0.95置信区间
344
    intervals = inference.calculate_confidence_intervals(candidates_dict, alpha)
345

346
    print("\n=== 估计结果 ===")
347
    for param, (lower, upper) in intervals.items():
348
        true_value = getattr(inference, f"{param.replace('_candidates', '')}_true", None)
349
        if true_value is not None:
350
            print(f"{param}: [{lower:.3f}, {upper:.3f}] (真实值: {true_value:.3f})")
351
        else:
352
            print(f"{param}: [{lower:.3f}, {upper:.3f}]")
353

354
    # 可视化结果
355
    print("\n生成可视化图表...")
356
    inference.plot_results(candidates_dict, intervals, Y_obs, W_obs)
357

358
    # 进行覆盖率模拟（样本量较小以加快演示速度）
359
    print("\n正在进行覆盖率模拟...")
360
    coverage_results, interval_lengths = inference.coverage_simulation(
361
        n_outer_trials=500,  # 减少试验次数以加快演示
362
        n_inner_samples=1000,  # 减少内部样本数
363
        alpha=0.05
364
    )
365

366
    print("\n=== 覆盖率结果 ===")
367
    for param, coverage in coverage_results.items():
368
        print(f"{param}的覆盖率: {coverage:.3f} (目标: {1-alpha:.3f})")
369

370
    # 绘制区间长度分布
371
    plt.figure(figsize=(12, 4))
372
    for i, (param, lengths) in enumerate(interval_lengths.items()):
373
        plt.subplot(1, 3, i+1)
374
        plt.hist(lengths, bins=30, alpha=0.7)
375
        plt.xlabel(f'{param}区间长度')
376
        plt.ylabel('频数')
377
        plt.title(f'{param}置信区间长度分布')
378
    plt.tight_layout()
379
    plt.show()
380
if __name__ == "__main__":
381
    main()

音乐

音乐