我有3个关于scipy.stats.rv_continuous子类的问题。我的目标是编码一个统计混合模型的截断正态分布,截断指数分布和2个均匀分布。
( 1)为什么通过mm_model.rvs(大小= 1000)绘制随机变量的速度非常慢?我在纪录片中读到了一些关于表演问题的文章,但我真的很惊讶。
2)通过mm_model.rvs(大小=1,000)绘制随机变量后,得到IntegrationWarning?
IntegrationWarning:最大的细分数(50)已经实现。如果增加限制没有得到改善,则建议对被积体进行分析,以确定困难。如果能够确定局部困难的位置(奇异性、不连续性),则很可能通过分割区间并调用子区间上的积分器来获得该位置。也许应该使用一个特殊用途的集成商。warnings.warn(msg,IntegrationWarning)
3)我在纪录片中看到,我可以通过“形状”参数将参数传递给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()发布于 2017-03-02 18:23:42
(1)和(2)可能是相关的。你是要求枕木根据你提供的密度来产生随机样本。
我真的不知道map做什么,但我怀疑它集成了密度("pdf"),得到概率函数("cdf"),然后反演它来映射均匀样本到您的分布。这在数字上是昂贵的,就像您所经历的那样容易出错。
为了加快速度,您可以通过直接实现_rvs来帮助参与。只需画一个均匀的,以决定您的混合物选择哪个子模型,然后调用所选子模型的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这里是枕骨医生要说的话。
关于形状的注意事项:子类不需要显式指定它们。在这种情况下,形状将自动从被覆盖的方法(pdf,cdf等)的签名中推导出来。如果出于某种原因,您希望避免依赖内省,则可以显式地将形状指定为实例构造函数的参数。
所以通常的方法是在你的方法中加入一些参数。
https://stackoverflow.com/questions/42552117
复制相似问题