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一些提示：

要批量选择子模型。由于您的子模型的数量很小，所以对此应该足够好。如果可能的话，对于子模型使用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 

Re（3）这是scipy文档必须说的。

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

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