Быстрая выборка из пользовательского дистрибутива

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