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

Question

在这篇文章之后，我尝试通过创建类来创建 logit 正态分布LogitNormal：

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但这是一个完全不同的故事。

Answer 1

前言

我本周遇到了这个问题，我发现的唯一相关问题是这篇文章。我和OP有几乎相同的要求：

• 具有Logit 正态分布的随机变量。

但我还需要：

• 也能够进行统计检验；
• 同时符合scipy随机变量接口。

正如所@Jacques Gaudin指出的，for 的接口rv_continous（详细信息请参阅分发架构）在继承此类时不确保后续的 forloc和参数。scale这在某种程度上是误导性的，也是不幸的。

实现该__init__方法当然允许创建丢失的绑定，但代价是：它破坏了当前用于实现随机变量的模式（请参阅 lognormalscipy的实现示例）。

因此，我花了一些时间深入研究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（这个样本很正常）如下所示：