python - 从用户定义的分布中快速采样

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