连续分布的 scipy.stats 属性“熵”不能手动工作

Posted 2023-03-12

【中文标题】连续分布的 scipy.stats 属性“熵”不能手动工作

【英文标题】：scipy.stats attribute `entropy` for continuous distributions doesn't work manually

【发布时间】：2021-04-14 20:28:25

【问题描述】：

scipy.stats 中的每个连续分布都带有一个计算其微分熵的属性：.entropy。与正态分布 (norm) 和其他具有熵的封闭形式解的其他分布不同，其他分布必须依赖于数值积分。

试图找出在这些情况下.entropy 属性调用的是哪个函数，我在scipy.stats._distn_infrastructure.py 中找到了一个名为_entropy 的函数，它与integrate.quad(pdf) 一起使用（数值积分）。

但是当我尝试比较这两种方法时（属性.entropy 与函数_entropy 的数值积分），函数给出错误：

AttributeError: 'rv_frozen' object has no attribute '_pdf'

为什么分布的属性.entropy计算正常，而函数_entropy却报错？

import numpy as np
from scipy import integrate 
from scipy.stats import norm, johnsonsu
from scipy.special import entr

def _entropy(self, *args): #from _distn_infrastructure.py
    def integ(x):
        val = self._pdf(x, *args)
        return entr(val)

    # upper limit is often inf, so suppress warnings when integrating
    # _a, _b = self._get_support(*args)
    _a, _b = -np.inf, np.inf   
    with np.errstate(over='ignore'):
        h = integrate.quad(integ, _a, _b)[0]

    if not np.isnan(h):
        return h
    else:
        # try with different limits if integration problems
        low, upp = self.ppf([1e-10, 1. - 1e-10], *args)
        if np.isinf(_b):
            upper = upp
        else:
            upper = _b
        if np.isinf(_a):
            lower = low
        else:
            lower = _a
    return integrate.quad(integ, lower, upper)[0]

使用该属性可以正常工作：

print(johnsonsu(a=2.55,b=2.55).entropy())

返回0.9503703091220894

但函数没有：

print(_entropy(johnsonsu(a=2.55,b=2.55)))

返回错误 AttributeError: 'rv_frozen' object has no attribute '_pdf'，即使johnsonsu does have this attribute：

def _pdf(self, x, a, b):
    # johnsonsu.pdf(x, a, b) = b / sqrt(x**2 + 1) *
    #                          phi(a + b * log(x + sqrt(x**2 + 1)))
    x2 = x*x
    trm = _norm_pdf(a + b * np.log(x + np.sqrt(x2+1)))
    return b*1.0/np.sqrt(x2+1.0)*trm

在johnsonsu的情况下，属性.entropy调用的是哪个函数？

【参考方案1】：

如果您使用的是冻结发行版，则需要johnsonsu(a=2.55,b=2.55).entropy()，否则需要johnsonsu.entropy(a=2.55,b=2.55)。

您问题的为什么部分基本上是 _entropy 中的前导下划线表示“实现细节，不要直接调用”。一个更长的答案是，冻结的发行版包装了一个发行版实例（self.dist），并将对 _pdf、_pmf 等的调用委托给它。

编辑：执行johnsonsu(a=2.55,b=2.55) 创建一个冻结的发行版rv_frozen。除非您想多次重用实例，否则不要这样做：只需将 a,b 形状参数作为参数提供给熵函数。

【讨论】：

那么我怎样才能让手动_entropy 函数工作呢？这真的是.entropy 属性所调用的吗？我不知道我是否使用了冻结发行版，因为我不知道它们是什么。我以为我在做的只是估计johnsonsu。您编写的两个代码示例对冻结和非冻结之间的区别没有帮助，因为它们彼此都可以正常工作