通过子类化scipy.stats.rv_continuous创建偏态正态分布的问题

Question

更新：我找到了分布的问题。现在大部分都能正常工作了，除了当形状参数为负数的时候。虽然PDF应该能处理负的形状值，但在我子类化的分布中却不行。

---

我正在尝试用scipy的stats库创建一个偏态正态分布。目前我只需要PDF。

我对rv_continuous进行了子类化，但当我使用skew_norm.pdf(x, shape)时，得到的是一堆NaN。

这是我的类：

class skew_norm_gen(rv_continuous):
    def _pdf(self, x, s):
        return 2 * norm.pdf(x) * norm.cdf(x * s)

skew_norm = skew_norm_gen(name='skew_norm', shapes='s')

我尝试在类外直接计算PDF，这样是可以的。

另外，如果我加上*args*，我能像对待正态分布PDF那样传递位置和尺度吗？比如norm.pdf(x, loc=mu, scale=std)：

class skew_norm_gen(rv_continuous):
    def _pdf(self, x, s, *args):
        return 2 * norm.pdf(x, *args) * norm.cdf(x * s, *args)

skew_norm = skew_norm_gen(name='skew_norm', shapes='s')

谢谢。

---

更新：

我还尝试了一个简单的例子，感谢CT Zhu的建议。下面的代码有时候会输出NaN数组，有时候又会输出一组值。

In [26]:
import scipy.stats as ss

class skew_norm_gen(ss.rv_continuous):
    def _pdf(self, x, s):
        return 2 * ss.norm.pdf(x) * ss.norm.cdf(x * s)
skew_norm = skew_norm_gen(name='skew_norm', shapes='s')

In [27]:
data = ss.norm.rvs(0, size=100)
s = ss.skew(data)
skew_norm.pdf(data, s)

Out[28]:
array([ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
        nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
        nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
        nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
        nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
        nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
        nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
        nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
        nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
        nan])

---

更新 2：

如果形状参数小于0，PDF会输出NaN。

我可以直接计算skewnorm的PDF，结果是正常的。如果我尝试使用子类化的PDF，它就返回NaNs。

Answer 1

无法重现这个错误，看看：

In [15]:
import scipy.stats as ss
class skew_norm_gen(ss.rv_continuous):
    def _pdf(self, x, s):
        return 2 * ss.norm.pdf(x) * ss.norm.cdf(x * s)
skew_norm = skew_norm_gen(name='skew_norm', shapes='s')

In [17]:
skew_norm.pdf(3, 4)
Out[17]:
0.0088636968238760151

是的，你可以传递额外的 *args：

In [18]:

class skew_norm_gen(ss.rv_continuous):
    def _pdf(self, x, s, *args):
        return 2 * ss.norm.pdf(x, *args) * ss.norm.cdf(x * s, *args)
skew_norm = skew_norm_gen(name='skew_norm', shapes='s')

In [20]:
skew_norm.pdf(3, 4, loc=0.5, scale=3)
Out[20]:
0.18786061213807126

In [21]:
skew_norm.pdf(3, s=4, loc=0.5, scale=3)
Out[21]:
0.18786061213807126
In [22]:

skew_norm.pdf(3, s=4, loc=0, scale=1)
Out[22]:
0.0088636968238760151
In [28]:
plt.plot(np.linspace(-5, 5), skew_norm.pdf(np.linspace(-5,5),4), label='Skewed')
plt.plot(np.linspace(-5, 5), ss.norm.pdf(np.linspace(-5,5)), label='Normal')
plt.legend()    
Out[28]:
[<matplotlib.lines.Line2D at 0x1092667d0>]

编辑：

在你的示例数据中，s 是负数，这导致生成的 pdf 只包含 nan，这是 rv_continuous 定义的默认 badvalue（我想这就是它的意思）。

问题的根源在于：有一个默认的 _argcheck() 方法，用来验证参数是否有效。默认情况下，它会检查所有参数是否大于 0。在这个例子中，结果并不是这样。

所以解决方案是重写默认的 _argcheck() 方法，方法如下：

class skew_norm_gen(ss.rv_continuous):
    def _argcheck(self, skew):
        return np.isfinite(skew) #I guess we can confine it to finite value
    def _pdf(self, x, skew):
        return 2 * ss.norm.pdf(x) * ss.norm.cdf(x * skew)  

然后它应该就能正常工作了。

（另外，我建议把额外的参数叫做 skew，这样更容易理解。's' 可能会让人联想到标准差等等。）