scipy.stats.rv_continuous 的子类化

Question

我有 3 个关于 scipy.stats.rv_continuous 的子类化的问题。 我的目标是编写截断正态分布、截断指数分布和 2 个均匀分布的统计混合模型。

1) 为什么通过 mm_model.rvs(size = 1000) 绘制随机变量如此缓慢？我在纪录片中读到了一些关于性能问题的内容，但我真的很惊讶。

2) 通过 mm_model.rvs(size = 1000) 绘制随机变量后，我得到了这个 IntegrationWarning？

IntegrationWarning：已达到最大细分数 (50)。 如果增加限制没有改善，建议分析 被积函数以确定困难。如果一个位置 可以确定局部难度（奇点、不连续性） 可能会从分割区间和调用积分器中获益 在子范围内。也许应该使用专用的积分器。 warnings.warn(msg, IntegrationWarning)

3）我在纪录片中读到我可以通过“shape”参数将参数传输到pdf。我试图调整我的 pdf 并设置形状参数，但它不起作用。谁能解释一下？

感谢您的帮助。

def trunc_norm(z,low_bound,up_bound,mu,sigma):
    a = (low_bound - mu) / sigma
    b = (up_bound - mu) / sigma
    return stats.truncnorm.pdf(z,a,b,loc=mu,scale=sigma)



def trunc_exp(z,up_bound,lam):
    return stats.truncexpon.pdf(z,b=up_bound*lam,scale=1/lam)



def uniform(z,a,b):
    return stats.uniform.pdf(z,loc=a,scale=(b-a))


class Measure_mixture_model(stats.rv_continuous):

    def _pdf(self,z):

        z_true = 8
        z_min = 0
        z_max = 10
        p_hit = 0.7
        p_short = 0.1   
        p_max = 0.1
        p_rand = 0.1
        sig_hit = 1
        lambda_short = 0.5

        sum_p = p_hit + p_short + p_max + p_rand

        if sum_p < 0.99 or 1.01 < sum_p:
            misc.eprint("apriori probabilities p_hit, p_short, p_max, p_rand have to be 1!")
            return None

        pdf_hit = trunc_norm(z,z_min,z_max,z_true,sig_hit)
        pdf_short = trunc_exp(z,z_true,lambda_short)
        pdf_max = uniform(z,0.99*z_max,z_max)
        pdf_rand = uniform(z,z_min,z_max)

        pdf = p_hit * pdf_hit + p_short * pdf_short + p_max * pdf_max + p_rand * pdf_rand

        return pdf




#mm_model = Measure_mixture_model(shapes='z_true,z_min,z_max,p_hit,p_short,p_max,p_rand,sig_hit,lambda_short')
mm_model = Measure_mixture_model()


z = np.linspace(-1,11,1000)

samples = mm_model.pdf(z)
plt.plot(z,samples)
plt.show()


rand_samples = mm_model.rvs(size = 1000)


bins = np.linspace(-1, 11, 100)
plt.hist(rand_samples,bins)
plt.title("measure mixture model")
plt.xlabel("z: measurement")
plt.ylabel("p: relative frequency")
plt.show()

Answer 1

(1) 和 (2) 可能是相关的。您要求 scipy 仅根据您提供的密度生成随机样本。

我真的不知道 scipy 做了什么，但我怀疑它整合了密度（“pdf”）以获得概率函数（“cdf”），然后将其反转以将均匀样本映射到您的分布。这在数值上很昂贵，而且您容易出错。

为了加快速度，您可以通过直接实现 _rvs 来帮助 scipy。只需绘制一个制服来决定选择混合的哪个子模型，然后调用所选子模型的rvs。您可能需要的其他功能也类似。

这里有一些关于如何实现矢量化rvs 的提示：

批量选择子模型。由于您的子模型数量很少np.digitize 应该足以满足此要求。如果可能的话，为子模型使用rv_frozen 实例；它们非常方便，但我似乎记得您不能将所有可选参数都传递给它们，因此您可能必须单独处理。

self._bins = np.cumsum([p_hit, p_short, p_max])
self._bins /= self._bins[-1] + p_rand

submodel = np.digitize(uniform.rvs(size=size), self._bins)

result = np.empty(size)
for j, frozen in enumerate((frz_trunc_norm, frz_trunc_exp, frz_unif_1, frz_unif_2)):
    inds = np.where(submodel == j)
    result[inds] = frozen.rvs(size=inds.shape)
return result

关于 (3) 这是 scipy 文档要说的内容。

关于形状的说明：子类不需要明确指定它们。在这种情况下，将自动从重写方法（pdf、cdf 等）的签名中推断出形状。如果出于某种原因，您希望避免依赖自省，则可以显式指定形状作为实例构造函数的参数。

所以通常的方法是在你的方法中加入一些参数。