计算概率分布的微分熵和信息熵的JVHW算法实现

共7个文件

py：4个

ipynb：2个

mat：1个

版权申诉

python

开发语言

109 浏览量 2022-04-01 15:31:24 上传评论收藏 36KB ZIP 举报

在IT领域，尤其是在数据分析、机器学习以及信息理论中，概率分布的熵是一个极其重要的概念。熵是用来衡量一个随机变量不确定性或信息含量的度量。在这个主题中，我们关注的是微分熵和信息熵，以及如何使用JVHW（Jensen-Van Hoesen-Wallace）算法来计算它们。这篇文章将探讨这两个概念及其Python实现。我们要理解信息熵。信息熵，最初由Claude Shannon提出，是用于衡量随机变量平均信息量的函数。对于离散随机变量X，其概率质量函数为p(x)，信息熵H(X)定义为： \[ H(X) = -\sum_{x} p(x) \log_2 p(x) \] 这个公式告诉我们，熵越大，随机变量的不确定性越高。微分熵则是信息熵在连续随机变量上的扩展。对于连续随机变量Y，其概率密度函数为f(y)，微分熵H(Y)定义为： \[ H(Y) = -\int_{-\infty}^{\infty} f(y) \log_2 f(y) dy \] 微分熵同样反映了连续随机变量的不确定性，但其积分形式更适应于连续分布。 JVHW算法是由Jensen、Van Hoesen、Wallace三位学者提出的，主要用于估计概率分布的熵。该算法基于插值和拟合方法，通过有限的样本数据来估计整个分布的熵。它适用于处理各种类型的数据，包括非均匀分布和具有复杂结构的数据。在Python中实现JVHW算法，我们可以利用科学计算库如NumPy、SciPy以及Pandas等。我们需要收集样本数据，然后对数据进行预处理，如归一化或者标准化。接着，可以使用插值技术（如Scipy的interpolate模块）来估计概率密度函数。然后，根据概率密度函数，计算微分熵。以下是基本步骤： 1. 导入所需库： ```python import numpy as np from scipy.interpolate import interp1d from scipy.stats import entropy import pandas as pd ``` 2. 收集并预处理数据： ```python data = np.array([...]) # 假设你的数据在这里 data = (data - data.min()) / (data.max() - data.min()) # 数据归一化 ``` 3. 插值得到概率密度函数： ```python f = interp1d(data, np.ones(len(data)), kind='cubic') # 使用立方插值 ``` 4. 计算微分熵： ```python pdf = f(data) / np.sum(f(data)) # 计算概率密度函数 differential_entropy = entropy(pdf, base=2) # 使用Scipy的entropy函数计算微分熵 ``` 以上就是一个基本的Python实现过程。在实际应用中，可能需要调整插值方法、考虑数据分布特性等因素以提高估计精度。此外，还可以使用其他方法，如KDE（Kernel Density Estimation）来估计概率密度函数。总结来说， JVHW算法为我们在Python中计算概率分布的微分熵和信息熵提供了一种有效的方法。通过理解和应用这些概念，我们可以更好地理解和分析数据的不确定性和复杂性，这对于数据驱动的决策和模型建立至关重要。

资源推荐

资源详情

资源评论

收起资源包目录

计算概率分布的微分熵和信息熵的JVHW算法实现_python.zip （7个子文件）

Python

MI_estimation_demo.ipynb 22KB

test_MI.py 2KB

est_MI.py 3KB

est_entro.py 5KB

test_entro.py 2KB

Entropy_estimation_demo.ipynb 22KB

poly_coeff_entro.mat 3KB

import numpy as np import scipy.io as sio poly_entro = None def est_entro_JVHW(samp): """Proposed JVHW estimate of Shannon entropy (in bits) of the input sample This function returns a scalar JVHW estimate of the entropy of samp when samp is a vector, or returns a row vector containing the JVHW estimate of each column of samp when samp is a matrix. Input: ----- samp: a vector or matrix which can only contain integers. The input data type can be any interger classes such as uint8/int8/ uint16/int16/uint32/int32/uint64/int64, or floating-point such as single/double. Output: ----- est: the entropy (in bits) of the input vector or that of each column of the input matrix. The output data type is double. """ samp = formalize_sample(samp) [n, wid] = samp.shape n = float(n) # The order of polynomial is no more than 22 because otherwise floating-point error occurs order = min(4 + int(np.ceil(1.2 * np.log(n))), 22) global poly_entro if poly_entro is None: poly_entro = sio.loadmat('poly_coeff_entro.mat')['poly_entro'] coeff = poly_entro[order-1, 0][0] f = fingerprint(samp) prob = np.arange(1, f.shape[0] + 1) / n # Piecewise linear/quadratic fit of c_1 V1 = np.array([0.3303, 0.4679]) V2 = np.array([-0.530556484842359, 1.09787328176926, 0.184831781602259]) f1nonzero = f[0] > 0 c_1 = np.zeros(wid) with np.errstate(divide='ignore', invalid='ignore'): if n >= order and f1nonzero.any(): if n < 200: c_1[f1nonzero] = np.polyval(V1, np.log(n / f[0, f1nonzero])) else: n2f1_small = f1nonzero & (np.log(n / f[0]) <= 1.5) n2f1_large = f1nonzero & (np.log(n / f[0]) > 1.5) c_1[n2f1_small] = np.polyval(V2, np.log(n / f[0, n2f1_small])) c_1[n2f1_large] = np.polyval(V1, np.log(n / f[0, n2f1_large])) # make sure nonzero threshold is higher than 1/n c_1[f1nonzero] = np.maximum(c_1[f1nonzero], 1 / (1.9 * np.log(n))) prob_mat = entro_mat(prob, n, coeff, c_1) return np.sum(f * prob_mat, axis=0) / np.log(2) def entro_mat(x, n, g_coeff, c_1): # g_coeff = {g0, g1, g2, ..., g_K}, K: the order of best polynomial approximation, K = len(g_coeff) - 1 thres = 4 * c_1 * np.log(n) / n T, X = np.meshgrid(thres, x) ratio = np.minimum(np.maximum(2 * X / T - 1, 0), 1) q = np.arange(K).reshape((1, 1, K)) g = g_coeff.reshape((1, 1, K + 1)) MLE = - X * np.log(X) + 1 / (2 * n) polyApp = np.sum(np.concatenate((T[..., None], ((n * X)[..., None] - q) / (T[..., None] * (n - q))), axis=2).cumprod(axis=2) * g, axis=2) - X * np.log(T) polyfail = np.isnan(polyApp) | np.isinf(polyApp) polyApp[polyfail] = MLE[polyfail] output = ratio * MLE + (1 - ratio) * polyApp return np.maximum(output, 0) def est_entro_MLE(samp): """Maximum likelihood estimate of Shannon entropy (in bits) of the input sample This function returns a scalar MLE of the entropy of samp when samp is a vector, or returns a (row-) vector consisting of the MLE of the entropy of each column of samp when samp is a matrix. Input: ----- samp: a vector or matrix which can only contain integers. The input data type can be any interger classes such as uint8/int8/ uint16/int16/uint32/int32/uint64/int64, or floating-point such as single/double. Output: ----- est: the entropy (in bits) of the input vector or that of each column of the input matrix. The output data type is double. """ samp = formalize_sample(samp) [n, wid] = samp.shape n = float(n) f = fingerprint(samp) prob = np.arange(1, f.shape[0] + 1) / n prob_mat = - prob * np.log2(prob) return prob_mat.dot(f) def formalize_sample(samp): samp = np.array(samp) if np.any(samp != np.fix(samp)): raise ValueError('Input sample must only contain integers.') if samp.ndim == 1 or samp.ndim == 2 and samp.shape[0] == 1: samp = samp.reshape((samp.size, 1)) return samp def fingerprint(samp): """A memory-efficient algorithm for computing fingerprint when wid is large, e.g., wid = 100 """ wid = samp.shape[1] d = np.r_[ np.full((1, wid), True, dtype=bool), np.diff(np.sort(samp, axis=0), 1, 0) != 0, np.full((1, wid), True, dtype=bool) ] f_col = [] f_max = 0 for k in range(wid): a = np.diff(np.flatnonzero(d[:, k])) a_max = a.max() hist, _ = np.histogram(a, bins=a_max, range=(1, a_max + 1)) f_col.append(hist) if a_max > f_max: f_max = a_max return np.array([np.r_[col, [0] * (f_max - len(col))] for col in f_col]).T

评论收藏

内容反馈

版权申诉