使用自定义概率密度函数时,stats.rv_continuous很慢

3

我希望能够可视化两个由数据估算出来的概率密度函数之间的联合分布(都是通过KDE获得)。假设有一个KDE,我有一个排序为元组data的离散x,y数据。我需要生成具有此分布的随机变量,以执行概率积分变换(最终获得均匀分布)。我生成随机变量的方法如下:

import scipy.stats as st
from scipy import interpolate, integrate

pdf1 = interpolate.interp1d(data[0], data[1])

class pdf1_class(st.rv_continuous):
    def _pdf(self,x):
        return pdf1(x)

pdf1_rv = pdf1_class(a = data[0][0], b= data[0][-1], name = 'pdf1_class')

pdf1_samples = pdf1_rv.rvs(size=10000)

然而,这种方法非常缓慢。我还收到了以下警告:

IntegrationWarning: 达到了最大细分数(50)。 如果增加限制没有改善,则建议分析被积函数以确定困难之处。 如果可以确定局部困难的位置(奇点、不连续性),则我们可能会从将区间分割并在子范围上调用积分器中获益。 可能应该使用专用积分器。 warnings.warn(msg, IntegrationWarning)

IntegrationWarning: 检测到舍入误差,这阻止了达到所请求的容差。误差可能被低估。 warnings.warn(msg, IntegrationWarning)

有更好的方法生成随机变量吗?

根据文档,可以通过覆盖 _logpdf,_cdf,_logcdf,_ppf,_rvs,_isf,_sf,_logsf 来提高速度。要加快 rvs,覆盖 _ppf(百分点函数,cdf 的反函数)可能已足够。 - unutbu
我覆盖了 _cdf,但没有什么显著的区别。然后我也覆盖了 _ppf ,现在10000个样本是瞬间完成的。谢谢你。你想把这作为答案吗?我会接受它的。 - mch56
既然您已经付出了实现 _ppf 的努力,我认为您将是撰写出色答案的最佳人选。(回答自己的问题不仅可以接受,而且是鼓励的。) - unutbu
1个回答

3

根据@unutbu的建议,我实现了_cdf_ppf,这使得10000个样本的计算变得瞬间完成。为此,我在上述代码中添加了以下内容:

discrete_cdf1 = integrate.cumtrapz(y=data[1], x = data[0])
cdf1 = interpolate.interp1d(data[0][1:], discrete_cdf1)
ppf1 = interpolate.interp1d(discerete_cdf1, data[0][:-1])

我随后在pdf1_class中添加了以下两个方法:

def _cdf(self,x):
    return cdf1(x)

def _ppf(self,x):
    return ppf1(x)

2
似乎应该跳过第一个x,而不是最后一个:cdf1 = interpolate.interp1d(data[0][1:], discrete_cdf1) - Anatoly
好的,会修改。 - mch56
1
你也可以在 cumtrapz 中使用 initial = 0 作为关键字参数。这样就不需要将数据向量的大小减少1了。 - mch56
2
同时对于ppf,ppf1 = interpolate.interp1d(discerete_cdf1, data[0][1:]) - Anatoly

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接