numpy Python计算熵的问题

xzv2uavs  于 2023-05-17  发布在  Python
关注(0)|答案(1)|浏览(171)

我试图计算微分熵(来自信息论),但在Python中遇到了一些问题。我的尝试如下:
我有下面的微分熵函数:

import numpy as np
from scipy.stats import norm
from scipy import integrate

def diff_entropy(nu, constant):

  def pdf_gaus_mixture(input):
    return (1-nu)*norm.pdf(input, loc=0, scale=1) + nu*norm.pdf(input, loc=constant, scale=1)

  def func(input):
    return pdf_gaus_mixture(input) * np.log(1 / pdf_gaus_mixture(input))

  return integrate.quad(func, -np.inf, np.inf)[0]

我想计算如下:

nu=0.1
beta=0.01
delta=0.1
sigma=0.01
diff_entropy(nu, np.sqrt(1/((beta/delta)+(sigma**2))))

但是Python给了我以下错误:

<ipython-input-22-6267f1f9e56a>:7: RuntimeWarning: divide by zero encountered in double_scalars
  return pdf_gaus_mixture(input) * np.log(1 / pdf_gaus_mixture(input))
<ipython-input-22-6267f1f9e56a>:7: RuntimeWarning: invalid value encountered in double_scalars
  return pdf_gaus_mixture(input) * np.log(1 / pdf_gaus_mixture(input))
<ipython-input-22-6267f1f9e56a>:7: RuntimeWarning: overflow encountered in double_scalars
  return pdf_gaus_mixture(input) * np.log(1 / pdf_gaus_mixture(input))
<ipython-input-22-6267f1f9e56a>:9: IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  return integrate.quad(func, -np.inf, np.inf)[0]
nan

**问题:**我做错了什么?我怀疑这个问题是由于积分的端点是负无穷大和正无穷大。我可以把它改得很小,比如正负10,但我担心近似值的准确性会有所损失。有没有更聪明的方法来克服这一点?谢谢。

fhg3lkii

fhg3lkii1#

嵌套函数func0.0乘以np.inf的不同倍数,这是未定义的。我发现这样修改你的函数:

def diff_entropy(nu, constant):

  def pdf_gaus_mixture(input):

    return (1-nu)*norm.pdf(input, loc=0, scale=1) + nu*norm.pdf(input, loc=constant, scale=1)

  def func(input):
    expr1 = pdf_gaus_mixture(input)
    expr2 = np.log(1 / pdf_gaus_mixture(input))

    if expr1 == 0:
        print(input, expr1, expr2)
    return expr1 * expr2

  return integrate.quad(func, -np.inf, np.inf)[0]

从技术上讲,您可以尝试循环计算并增加积分的下限和上限,直到python将0乘以np.inf,也就是说,直到python无法给予您更准确的结果。我使用下面的代码来实现这一点。如果有用的话告诉我。

import numpy as np
from scipy.stats import norm
from scipy import integrate

def diff_entropy(nu, constant, lower_int_boundary, upper_int_boundary):

    def pdf_gaus_mixture(input):

        return (1-nu)*norm.pdf(input, loc=0, scale=1) + nu*norm.pdf(input, loc=constant, scale=1)

    def func(input):
        return pdf_gaus_mixture(input) * np.log(1 / pdf_gaus_mixture(input))

    return integrate.quad(func, lower_int_boundary, upper_int_boundary)[0]

nu=0.1
beta=0.01
delta=0.1
sigma=0.01
constant = np.sqrt(1/((beta/delta)+(sigma**2)))

lower_int_boundary = 0
upper_int_boundary = 0
step_size = 0.25
entropy_results = list()
boundaries = list()

while True:
    lower_int_boundary -= step_size
    upper_int_boundary += step_size

    entropy = diff_entropy(nu, constant, lower_int_boundary, upper_int_boundary)

    if np.isnan(entropy):
        break

    entropy_results.append(entropy)
    boundaries.append([lower_int_boundary, upper_int_boundary])

print(f"Most accurate entropy calculated: {entropy_results[-1]}")  # 1.6664093342815425
print(f"Boundaries used: {boundaries[-1]}")  # [-37.5, 37.5]

相关问题