如何在Python中创建一个logit-normal分布？

Question

按照这篇帖子，我尝试通过创建 LogitNormal 类来创建一个logit-normal分布：

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-normal分布的唯一原因是因为我找到的唯一实现该分布的方式是来自于PyMC3包和TensorFlow包，两者都相当沉重/过度，如果你只需要一个函数，那么它们就显得有些臃肿。我已经尝试过PyMC3，但似乎它对scipy的处理不好，它总是崩溃。但那是完全不同的故事了。

- mapf

---

为什么你要调用 LogitNormal().pdf？我猜你应该调用 LogitNormal()._pdf。因为 _pdf 是你在 LogitNormal 类中唯一方法的调用方式。 - user10325516

因为根据文档 https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.html 和我在问题中链接的帖子，重新定义 _pdf 或 _cdf 是如何对 rv_continuous 进行子类化的。 - mapf

根据源代码 rv_continuous.pdf，该方法使用 loc 和 scale 参数对 x 值进行转换。我猜你可以使用默认值得到相同的结果：sigma, mu = 1, 0。 - user10325516

我明白了。我有点困惑，这是什么意思？ - mapf

假设问题源于重复计算。rv_continuous.pdf计算概率密度函数。而scipy.stats.norm.pdf也计算概率密度函数。我想这就是我能告诉你的全部了，因为我不熟悉这些函数。祝你好运。 - user10325516

Answer 1

前言

本周我遇到了这个问题，唯一与之相关的问题是这篇帖子。我的需求几乎与OP相同：

• 有一个逻辑正态分布的随机变量。

但我还需要：

• 能够执行统计检验；
• 同时符合scipy随机变量接口的规范。

正如@Jacques Gaudin所指出的，rv_continous的接口（详见分布结构）并没有保证从该类继承时要遵循loc和scale参数，而这在某种程度上是误导性的和不幸的。

当然，通过实现__init__方法可以创建缺失的绑定，但这是一个折衷方案：它打破了scipy目前用于实现随机变量的模式（其中一个实现的例子为对数正态分布）。

因此，我花时间深入研究了scipy代码，并为这个分布创建了一个MCVE。虽然它不是完全完整的（主要缺失了矩重写），但它满足OP和我的需求，同时具有令人满意的精度和性能。

MCVE

这个随机变量的接口遵循规范的实现可以是：

class logitnorm_gen(stats.rv_continuous):

    def _argcheck(self, m, s):
        return (s > 0.) & (m > -np.inf)
    
    def _pdf(self, x, m, s):
        return stats.norm(loc=m, scale=s).pdf(special.logit(x))/(x*(1-x))
    
    def _cdf(self, x, m, s):
        return stats.norm(loc=m, scale=s).cdf(special.logit(x))
    
    def _rvs(self, m, s, size=None, random_state=None):
        return special.expit(m + s*random_state.standard_normal(size))
    
    def fit(self, data, **kwargs):
        return stats.norm.fit(special.logit(data), **kwargs)

logitnorm = logitnorm_gen(a=0.0, b=1.0, name="logitnorm")

这个实现可以释放大部分scipy随机变量的潜力。

N = 1000
law = logitnorm(0.24, 1.31)             # Defining a RV
sample = law.rvs(size=N)                # Sampling from RV
params = logitnorm.fit(sample)          # Infer parameters w/ MLE
check = stats.kstest(sample, law.cdf)   # Hypothesis testing
bins = np.arange(0.0, 1.1, 0.1)         # Bin boundaries
expected = np.diff(law.cdf(bins))       # Expected bin counts

由于它依赖于scipy正态分布，因此我们可以假设其底层函数具有与正态随机变量对象相同的准确性和性能。但是，当处理高度倾斜的分布时，特别是在支持边界处，它可能确实会受到浮点算术不精确的影响。

测试

为了检查它的表现如何，我们绘制一些感兴趣的分布并进行检查。 让我们创建一些固定装置:

def generate_fixtures(
    locs=[-2.0, -1.0, 0.0, 0.5, 1.0, 2.0],
    scales=[0.32, 0.56, 1.00, 1.78, 3.16],
    sizes=[100, 1000, 10000],
    seeds=[789, 123456, 999999]
):
    for (loc, scale, size, seed) in itertools.product(locs, scales, sizes, seeds):
        yield {"parameters": {"loc": loc, "scale": scale}, "size": size, "random_state": seed}

并且对相关的分布和样本进行检查：

eps = 1e-8
x = np.linspace(0. + eps, 1. - eps, 10000)

for fixture in generate_fixtures():
    
    # Reference:
    parameters = fixture.pop("parameters")
    normal = stats.norm(**parameters)
    sample = special.expit(normal.rvs(**fixture))
    
    # Logit Normal Law:
    law = logitnorm(m=parameters["loc"], s=parameters["scale"])
    check = law.rvs(**fixture)
    
    # Fit:
    p = logitnorm.fit(sample)
    trial = logitnorm(*p)
    resample = trial.rvs(**fixture)
    
    # Hypothetis Tests:
    ks = stats.kstest(check, trial.cdf)
    bins = np.histogram(resample)[1]
    obs = np.diff(trial.cdf(bins))*fixture["size"]
    ref = np.diff(law.cdf(bins))*fixture["size"]
    chi2 = stats.chisquare(obs, ref, ddof=2)

以下展示了使用n=1000, seed=789（此样本相当普通）进行的一些调整：

- jlandercy

Answer 2

如果您查看pdf方法的源代码，您会注意到在没有使用“scale”和“loc”关键字参数的情况下调用了“_pdf”。

   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()

- Jacques Gaudin

---

谢谢，这个有效！我知道这是一个不同类型的问题，但您是否也知道如果实现了反函数cdf / ppf会是什么样子呢？因为根据文档：“默认方法_rvs依赖于应用于均匀随机变量的cdf的反函数_ppf。为了有效地生成随机变量，需要覆盖默认_ppf（例如，如果可以明确表达反函数cdf）或在自定义_rvs方法中实施采样方法。” - mapf

我不了解反向 cdf / pdf，对统计学不是很了解。也许在 https://math.stackexchange.com/ 上询问会有帮助？ - Jacques Gaudin

好的，这里有一件奇怪的事情。我找到了如何做到这一点，并想相应地编辑您的答案，但界面不允许我这样做，声称代码格式不正确。然而，所谓的错误不在我的编辑中，而在您的引用中。这怎么可能呢？如果您直接打开编辑并尝试按原样保存，它将无法工作，除非您删除引用。 - mapf

奇怪，也许只需将其发布为您的答案。回答您的问题没有任何问题！ - Jacques Gaudin