如何以比这更快的速度在Python中采样非均质Poisson进程? [英] How to sample inhomogeneous Poisson processes in Python faster than this?

问题描述

我正在以不固定速率的毫秒级时间对Poisson进程进行采样.我通过检查每个大小增量的间隔中是否存在事件(基于该间隔中的平均速率)来离散化采样过程.由于我正在使用Python，因此它的运行速度比我希望的要慢一些.我当前使用的代码如下:

I'm sampling a Poisson process at a millisecond time scale where the rate is not fixed. I discretise the sampling process by checking in each interval of size delta whether there is an event there or not based on the average rate in that interval. Since I'm using Python it's running a bit slower than I would hope it to be. The code I'm currently using is the following:

import numpy
def generate_times(rate_function,max_t,delta):
    times = []
    for t in numpy.arange(delta,max_t,delta):
        avg_rate = (rate_function(t)+rate_function(t+delta))/2.0
        if numpy.random.binomial(1,1-math.exp(-avg_rate*delta/1000.0))>0:
            times.extend([t])
    return times

费率函数可以是任意的，我不是在给定费率函数的情况下寻找封闭形式的解决方案.

The rate function can be arbitrary, I'm not looking for a closed form solution given a rate function.

如果您想使用一些参数，可以尝试:

If you want some parameters to play with you can try:

max_t = 1000.0
delta = 0.1
high_rate = 100.0
low_rate = 0.0
phase_length = 25.0
rate_function = (lambda x: low_rate + (high_rate-low_rate)*0.5*(1+math.sin(2*math.pi*x/phase_length)))

推荐答案

以下是该版本的运行速度提高了约75倍，并实现了相同的功能:

Here's a version which runs about 75x faster and implements the same function:

def generate_times_opt(rate_function,max_t,delta):
    t = numpy.arange(delta,max_t, delta)
    avg_rate = (rate_function(t) + rate_function(t + delta)) / 2.0
    avg_prob = 1 - numpy.exp(-avg_rate * delta / 1000.0)
    rand_throws = numpy.random.uniform(size=t.shape[0])

    return t[avg_prob >= rand_throws].tolist()

输出:

11:59:07 [70]: %timeit generate_times(rate_function, max_t, delta)
10 loops, best of 3: 75 ms per loop

11:59:23 [71]: %timeit generate_times_opt(rate_function, max_t, delta)
1000 loops, best of 3: 1.15 ms per loop

---

边注:，但这不是模拟非均匀泊松过程的最佳方法.来自维基百科:

---

Sidenote: This is not the best way to simulate an Inhomogenous Poisson Process, though. From Wikipedia:

具有强度函数λ(t)的不均匀泊松过程可以通过从固定费率λ的均匀泊松过程中进行拒绝采样来模拟:选择一个足够大的λ，使得λ(t)=λp(t)并模拟一个费率参数为λ的泊松过程.在时间t处以概率p(t)接受来自Poisson模拟的事件.

An inhomogeneous Poisson process with intensity function λ(t) can be simulated by rejection sampling from a homogeneous Poisson process with fixed rate λ: choose a sufficiently large λ so that λ(t) = λ p(t) and simulate a Poisson process with rate parameter λ. Accept an event from the Poisson simulation at time t with probability p(t).

这篇关于如何以比这更快的速度在Python中采样非均质Poisson进程?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！