從用戶定義的分佈中快速採樣

Question

我運行的模擬需要我從概率分佈中繪製值。我這樣做如下：

import numpy as np 
import scipy.special as sp 
from scipy.stats import rv_continuous 


class epsilon_pdf(rv_continuous): 

    def _pdf(self, x, omega): 

     return np.exp(omega ** -1 * np.cos(x))/(2 * np.pi * 
                sp.iv(0, omega ** -1)) 


random_epsilon = epsilon_pdf(a=-np.pi, b=np.pi) 
n_trials = 1 # 10 ** 6 
goal_dict = {'omega': 2 ** -4, 'xi': 2 ** 0} 
for trial_num in xrange(n_trials): 
    # Choose m. 
    m = np.random.poisson(goal_dict['xi']) 
    # Draw m values for epsilon. 
    epsilon_values = random_epsilon.rvs(omega=goal_dict['omega'], size=m) 

（什麼是上面寫的是一個很小的玩具爲例）。

我遇到的一個主要問題是，調用random_epsilon.rvs是慢得令人難以置信 - 這麼慢當我將n_trials設置爲所需的10 ** 6時，某些值omega和xi會導致腳本花費377小時完成。

任何人都可以想到在我的概率分佈的Python代碼重新配置和我的採樣會更快嗎？ （也許有辦法做到這一點使用numpy的，這將是更快？）

（我不知道我的分佈是否已被命名標準之一。）

Answer 1

import numpy as np 
from scipy.stats import rv_continuous 
import scipy.special as sp 
import time 
from numpy import ma as ma 

class epsilon_pdf(rv_continuous): 
    def _pdf(self, x, omega): 
     return np.exp(omega ** -1 * np.cos(x))/(2 * np.pi * 
                sp.iv(0, omega ** -1)) 

random_epsilon = epsilon_pdf(a=-np.pi, b=np.pi) 
n_trials = 1000 # 10 ** 6 
goal_dict = {'omega': 2 ** -4, 'xi': 2 ** 0} 

random_epsilon_rvs = lambda x: random_epsilon.rvs(
    omega=goal_dict['omega'], 
    size = x 
    ) 

random_epsilon_rvs = np.vectorize(random_epsilon_rvs, otypes=[np.ndarray]) 

t0 = time.time() 
for trial_num in xrange(n_trials): 
    # Choose m. 
    m = np.random.poisson(goal_dict['xi']) 
    # Draw m values for epsilon. 
    epsilon_values = random_epsilon.rvs(omega=goal_dict['omega'], size=m) 

t1 = time.time() 
a = np.random.poisson(2**-4, size = (1, n_trials))[0] 
mask = a>0 
b = a[mask] 
c = random_epsilon_rvs(b) 

t2 = time.time() 

print t1-t0,t2-t1 
# 11.7730000019 vs 0.682999849319