Python scipy rv_continuous 实现的问题答案

【问题标题】：Issues with Python scipy rv_continuous implementationPython scipy rv_continuous 实现的问题
【发布时间】：2020-10-12 10:09:55
【问题描述】：

我正在尝试使用自定义分布创建 rv_continuous 的子类，我可以通过多个函数计算其 pdf。

这是我到目前为止所做的

import numpy as np
from scipy.stats import rv_continuous

辅助功能

def func1(xx, a_, b_, rho, m, sigma):
    return a_ + b_*(rho*(xx-m) + np.sqrt((xx-m)*(xx-m) + sigma*sigma))

def func2(xx, a_, b_, rho, m, sigma):
    sig2 = sigma*sigma
    return b_*(rho*np.sqrt((xx-m)*(xx-m)+sig2)+xx-m)/(np.sqrt((xx-m)*(xx-m)+sig2))

def func3(xx, a_, b_, rho, m, sigma):
    sig2 = sigma*sigma
    return b_*sig2/(np.sqrt((xx-m)*(xx-m)+sig2)*((xx-m)*(xx-m)+sig2))

def func4(xx, a_, b_, rho, m, sigma):
    w = func1(xx, a_, b_, rho, m, sigma)
    w1 = func2(xx, a_, b_, rho, m, sigma)
    w2 = func3(xx, a_, b_, rho, m, sigma)
    return (1.-0.5*xx*w1/w)*(1.0-0.5*xx*w1/w) - 0.25*w1*w1*(0.25 + 1./w) + 0.5*w2

def func5(xx, a_, b_, rho, m, sigma):
    vsqrt = np.sqrt(func1(xx, a_, b_, rho, m, sigma))
    return -xx/vsqrt - 0.5*vsqrt

最终的密度函数

def density(xx, a_, b_, rho, m, sigma):
    dm = func5(xx, a_, b_, rho, m, sigma)
    return func4(xx, a_, b_, rho, m, sigma)*np.exp(-0.5*dm*dm)/np.sqrt(2.*np.pi*func1(xx, a_, b_, rho, m, sigma))

一组参数

Params = 1.0073, 0.3401026, -0.8, 0.000830, 0.5109564

从函数中检查 pdf

xmin, xmax, nbPoints = -10., 10., 2000
x_real = np.linspace(xmin, xmax, nbPoints)

den_from_func = density(x_real, *Params)

现在构建我的分发类

class density_gen(rv_continuous):
    def _pdf(self, x, a_hat, b_hat, rho, m, sigma):
        return density(x, a_hat, b_hat, rho, m, sigma)

实例化

my_density = density_gen(name='density_gen')

my_density.a, my_density.b, my_density.numargs

正如我指定的 _pdf 我应该有一个工作分发实例

这行得通

pdf = my_density._pdf(x_real, *Params)

cdf 也可以工作，虽然它非常慢

cdf = my_density._cdf(x_real, *Params)
my_density._cdf(0.1, *Params)

但是对于所有其他方法，例如，我得到了 nans

my_density.mean(*Params)    
my_density.ppf(0.01, *Params)

我在这里做错了什么？

【问题讨论】：

你能用rvs生成随机数吗？如果不是，那么这可能是 SciPy 的一个限制，因为我可以在我自己的 PDF 采样器实现中使用您的 PDF 生成随机数。您提供的 PDF 是您的应用程序将使用的唯一 PDF 吗？
不，我不能。尝试了 my_density.rvs(100000, *Params) 并得到了 ValueError: Domain error in arguments。使用 my_density._rvs(100000, *Params) 给出 AttributeError: 'density_gen' 对象没有属性 '_size'。我也可以使用我自己的采样器从 pdf 生成样本，然后得到最适合 Scipy 的分布，但我真的很想有一个工作的 rv_continuous 实例。我每天（可能在日内）生成许多参数集，并且在许多情况下，pdf 可以有很大的不同。界限也可以改变，尽管程度较小。所以速度在这里是一个问题。谢谢！

标签： python random distribution scipy.stats

【解决方案1】：

看来您需要将_argcheck 方法添加到density_gen，因为您的分发使用自定义参数：

class density_gen(rv_continuous):

    def _argcheck(self, *Params):
        return True

    def _pdf(self, x, a_hat, b_hat, rho, m, sigma):
        return density(x, a_hat, b_hat, rho, m, sigma)

my_density = density_gen(name='density_gen')
pdf = my_density._pdf(x_real, *Params)
print(my_density.rvs(size=5, *Params))
print(my_density.mean(*Params))  
print(my_density.ppf(0.01, *Params))

但是rvs、mean等等会很慢，大概是因为该方法每次需要生成随机数或计算统计量时都需要对PDF进行积分。如果速度非常重要，您将因此需要向density_gen 添加一个使用自己的采样器的_rvs 方法。这方面的一个例子是我自己的DensityInversionSampler，当仅给出 PDF 和采样域时，它通过数值反转生成随机数。

【讨论】：

明白，如果我想在生产中使用它肯定需要添加更多的_methods。