如何在 Python 中创建 logit 正态分布？

Question

在this post 之后，我尝试通过创建LogitNormal 类来创建logit 正态分布：

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import logit
from scipy.stats import norm, rv_continuous


class LogitNormal(rv_continuous):
    def _pdf(self, x, **kwargs):
        return norm.pdf(logit(x), **kwargs)/(x*(1-x))


class OtherLogitNormal:
    def pdf(self, x, **kwargs):
        return norm.pdf(logit(x), **kwargs)/(x*(1-x))


fig, ax = plt.subplots()
values = np.linspace(10e-10, 1-10e-10, 1000)
sigma, mu = 1.78, 0
ax.plot(
    values, LogitNormal().pdf(values, loc=mu, scale=sigma), label='subclassed'
)
ax.plot(
    values, OtherLogitNormal().pdf(values, loc=mu, scale=sigma),
    label='not subclassed'
)
ax.legend()
fig.show()

但是，LogitNormal 类不会产生预期的结果。当我不继承 rv_continuous 时，它可以工作。这是为什么？我需要子类化才能工作，因为我还需要它附带的其他方法，例如rvs。

顺便说一句，我在 Python 中创建自己的 logit 正态分布的唯一原因是因为我能找到的该分布的唯一实现来自 PyMC3 package 和 TensorFlow package，两者都很重/ 如果你只需要它们来完成那一项功能，那就大材小用了。我已经尝试过PyMC3，但显然它不适用于scipy，我想，它总是让我崩溃。但这是一个完全不同的故事。

Answer 1

如果您查看source code of the pdf method，您会注意到在调用_pdf 时没有scale 和loc 关键字参数。

   if np.any(cond):
        goodargs = argsreduce(cond, *((x,)+args+(scale,)))
        scale, goodargs = goodargs[-1], goodargs[:-1]
        place(output, cond, self._pdf(*goodargs) / scale)

这会导致覆盖 _pdf 方法中的 kwargs 始终是一个空字典。

如果您仔细查看代码，您还会注意到缩放和位置由pdf 处理，而不是_pdf。

在您的情况下，_pdf 方法调用norm.pdf，因此loc 和scale 参数必须以某种方式在LogitNormal._pdf 中可用。

例如，您可以在创建LogitNormal 的实例时传递scale 和loc 并将值存储为类属性：

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import logit
from scipy.stats import norm, rv_continuous


class LogitNormal(rv_continuous):
    def __init__(self, scale=1, loc=0):
        super().__init__(self)
        self.scale = scale
        self.loc = loc

    def _pdf(self, x):
        return norm.pdf(logit(x), loc=self.loc, scale=self.scale)/(x*(1-x))


fig, ax = plt.subplots()
values = np.linspace(10e-10, 1-10e-10, 1000)
sigma, mu = 1.78, 0
ax.plot(
    values, LogitNormal(scale=sigma, loc=mu).pdf(values), label='subclassed'
)
ax.legend()
fig.show()