如何用Python正确模拟幂律型非齐次泊松过程(NHPP)？

你原代码使用稀疏（接受-拒绝）法模拟NHPP，存在4个致命问题，同时你对幂律NHPP的理论规律记忆存在偏差：

• 稀疏法要求取速率函数在整个模拟区间的最大值λ作为齐次泊松的模拟速率，你代码里把0到T的速率积分值（区间总期望事件数）当成了λ，完全不符合算法要求
• 未做时间截断：你设定的模拟上限T=500，但输出里最后一个事件时间到了2077，超出模拟区间的结果完全无效
• 接受判定逻辑错误：生成齐次泊松间隔的随机数和接受判定用的随机数必须是两个独立的均匀分布样本，你复用了同一个u，判定逻辑完全失效
• 规律认知偏差：幂律NHPP的间隔变化规律和你描述的正好相反：
  β>1：λ(t)随t递增，事件越来越密，相邻间隔递减
  β=1：λ(t)为常数，即齐次泊松过程，间隔平均长度恒定
  β<1：λ(t)随t递减，事件越来越疏，相邻间隔递增
  你测试β=0.5时期待看到间隔递减，本身就是错误预期。

幂律速率的NHPP（又称威布尔过程）的累积速率函数存在解析解，用逆变换法模拟比稀疏法效率更高，无接受损耗，结果完全精确：
累积速率函数 $\Lambda(t) = \int_0^t \lambda\beta s^{\beta-1}ds = \lambda t^\beta$，其逆函数为 $\Lambda^{-1}(x) = (x/\lambda)^{{1/\beta}$。根据NHPP模拟的逆变换原理：若$E_i$是速率为1的齐次泊松过程的第i个事件时间，则目标NHPP的第i个事件时间为$T_i=\Lambda}{-1}(E_i)$。

完整可运行代码如下：

import numpy as np
import matplotlib.pyplot as plt

def nhpp_powerlaw(beta, lam, T):
    """
    模拟幂律速率非齐次泊松过程
    参数：
        beta: 幂律指数，取值>0
        lam: 速率系数，取值>0
        T: 模拟时间上限，取值>0
    返回：
        event_times: 所有落在[0,T]区间内的事件时间，升序排列
        cum_counts: 对应事件时间的累计事件数，用于绘制N(t)曲线
    """
    event_times = []
    cum_e = 0  # 单位齐次泊松过程的累积时间
    while True:
        # 生成单位指数分布间隔
        cum_e += np.random.exponential(scale=1)
        t_i = (cum_e / lam) ** (1 / beta)
        if t_i > T:
            break
        event_times.append(t_i)
    cum_counts = np.arange(1, len(event_times)+1)
    return event_times, cum_counts

# 测试绘图，对应三个beta参数
if __name__ == "__main__":
    lam = 0.35
    T = 500
    param_set = [
        (0.5, "black", "β=0.5"),
        (1, "blue", "β=1"),
        (1.5, "red", "β=1.5")
    ]
    plt.figure(figsize=(10,6))
    for beta, color, label in param_set:
        t, n = nhpp_powerlaw(beta, lam, T)
        plt.step(t, n, color=color, label=label, where="post")
    plt.xlabel("t")
    plt.ylabel("N(t)")
    plt.legend()
    plt.grid(alpha=0.3)
    plt.show()

运行代码后输出的N(t)曲线完全匹配理论特征：

• β=0.5（黑线）：前期斜率高（事件密集、间隔短），后期斜率持续降低（事件稀疏、间隔变长），曲线上凸
• β=1（蓝线）：斜率恒定，为直线，和齐次泊松过程特征一致
• β=1.5（红线）：前期斜率低（事件稀疏、间隔长），后期斜率持续升高（事件密集、间隔变短），曲线下凸

内容的提问来源于stack exchange，提问作者Ruben Martins