p = np.linspace(1e-6, 1-1e-6, 100)
entropy = (p - 1) * np.log(1 - p) - p * np.log(p)
plt.figure(figsize=(4,4))
plt.plot(p, entropy)
plt.xlabel('p')
plt.ylabel('Shannon entropy om nats')
plt.show()
import math
def H(sentence):
"""
最优编码长度
"""
entropy = 0
# 这里有 256 个可能的 ASCII 符号
for character_i in range(256):
Px = sentence.count(chr(character_i)) / len(sentence)
if Px > 0:
entropy += -Px * math.log(Px, 2) # log 以 2 为底
return entropy
import random
import math
# 只用 64 个字符
simple_message = "".join([chr(random.randint(0, 64)) for i in range(500)])
print(simple_message)
H(simple_message)
'''
,?&.%>?'*?',?%@1
%.'4+9041=>?- ,@5.> 6(:4:;5#7:?3*@@#-'"$! 8,&.%/!'-9?@ 8%$ 8):+1%$9<
8.!=/)9217990+@<
1:=<1 =1:5=&#=1*6#6-57763# 4=*?*9@!<317)%*$,;1:">
5.933876645029648
'''
# KL 定义
from scipy.stats import entropy # 内置 KL
def kl(p, q):
"""
D(P || Q)
"""
p = np.asarray(p, dtype=np.float)
q = np.asarray(q, dtype=np.float)
return np.sum(np.where(p != 0, p * np.log(p / q), 0))
# 测试
p = [0.1, 0.9]
q = [0.1, 0.9]
print(entropy(p, q) == kl(p, q))
'''
True
'''
from scipy.stats import norm
# D(P || Q) 与 D(Q || P) 比较
x = np.linspace(1, 8, 500)
y1 = norm.pdf(x, 3, 0.5)
y2 = norm.pdf(x, 6, 0.5)
p = y1 + y2 # 构建 p(x)
KL_pq, KL_qp = [], []
q_list = []
for mu in np.linspace(0, 10, 500):
for sigma in np.linspace(0.1, 5, 50): # 寻找最优 q(x)
q = norm.pdf(x, mu, sigma)
q_list.append(q)
KL_pq.append(entropy(p, q))
KL_qp.append(entropy(q, p))
KL_pq_min = np.argmin(KL_pq)
KL_qp_min = np.argmin(KL_qp)
fig, axes = plt.subplots(1, 2, figsize=(10,3))
axes[0].set_ylim(0, 0.8)
axes[0].plot(x, p/2, 'b', label='$p(x)$')
axes[0].plot(x, q_list[KL_pq_min], 'g--', label='$q^*(x)$')
axes[0].set_xlabel('$x$')
axes[0].set_ylabel('$p(x)$')
axes[0].set_title('$q^*={arg\min}_ q D_ {KL}(p||q)$')
axes[1].set_ylim(0, 0.8)
axes[1].plot(x, p/2, 'b', label='$p(x)$')
axes[1].plot(x, q_list[KL_qp_min], 'g--', label='$q^*(x)$')
axes[1].set_xlabel('$x$')
axes[1].set_ylabel('$p(x)$')
axes[1].set_title('$q^*={arg\min}_ q D_ {KL}(q||p)$')
for ax in axes:
ax.legend(loc='upper right')