2015-09-22 52 views
3

我正在採樣一個毫秒時間尺度的泊松過程,其中速率不固定。我通過檢查尺寸增量的每個間隔來檢查抽樣過程是否存在或不存在基於該間隔內的平均速率的事件。由於我使用Python,它的運行速度比我希望的要慢。目前我正在使用的代碼如下:如何比Python更快地採樣Python中的不均勻泊松過程?

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 

率函數可以是任意的,我不找給定率函數封閉形式的解決方案。

如果你想要一些參數和你一起玩可以試試:

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))) 
+0

你在問怎麼優化''''generate_times'''嗎?它可以返回一個ndarray還是它必須是一個列表? – wwii

+0

它不需要完全優化,運行時的任何因子大於5都會產生巨大的差異。輸出必須是一個列表(如果速度更快......:P,則可以將ndarray轉換爲列表)。 – Haffi112

+0

'''generate_times'''似乎不完整,它從不向'''tunes'''增加任何東西。 – wwii

回答

4

這裏是它運行約75倍速度更快,實現相同功能的版本:

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 

旁註:雖然這不是模擬非均勻泊松過程的最佳方法。從Wikipedia

與強度函數λ(t)可以通過拒絕從與固定率λ一個齊次泊松過程採樣模擬不均勻的泊松過程:選擇一個足夠大的λ以便爲λ(t)=λp (t)並用速率參數λ模擬泊松過程。在概率p(t)的時間t接受來自泊松模擬的事件。

+0

太棒了,非常感謝,來自維基百科頁面的另一種選擇也很有意義。對於速率函數有一個最大值而其他值接近於0是弱的。 – Haffi112