练习
class MaxEntropy:
def init(self, EPS=0.005):
self._samples = []
self._Y = set() # 标签集合，相当去去重后的y
self._numXY = {} # key为(x,y)，value为出现次数
self._N = 0 # 样本数
self.Ep = [] # 样本分布的特征期望值
self._xyID = {} # key记录(x,y),value记录id号
self._n = 0 # 特征键值(x,y)的个数
self._C = 0 # 最大特征数
self._IDxy = {} # key为(x,y)，value为对应的id号
self._w = []
self._EPS = EPS # 收敛条件
self._lastw = [] # 上一次w参数值
def loadData(self, dataset): self._samples = deepcopy(dataset) for items in self._samples: y = items[0] X = items[1:] self._Y.add(y) # 集合中y若已存在则会自动忽略 for x in X: if (x, y) in self._numXY: self._numXY[(x, y)] += 1 else: self._numXY[(x, y)] = 1 self._N = len(self._samples) self._n = len(self._numXY) self._C = max([len(sample) - 1 for sample in self._samples]) self._w = [0] * self._n self._lastw = self._w[:] self._Ep_ = [0] * self._n for i, xy in enumerate(self._numXY): # 计算特征函数fi关于经验分布的期望 self._Ep_[i] = self._numXY[xy] / self._N self._xyID[xy] = i self._IDxy[i] = xy def _Zx(self, X): # 计算每个Z(x)值 zx = 0 for y in self._Y: ss = 0 for x in X: if (x, y) in self._numXY: ss += self._w[self._xyID[(x, y)]] zx += math.exp(ss) return zx def _model_pyx(self, y, X): # 计算每个P(y|x) zx = self._Zx(X) ss = 0 for x in X: if (x, y) in self._numXY: ss += self._w[self._xyID[(x, y)]] pyx = math.exp(ss) / zx return pyx def _model_ep(self, index): # 计算特征函数fi关于模型的期望 x, y = self._IDxy[index] ep = 0 for sample in self._samples: if x not in sample: continue pyx = self._model_pyx(y, sample) ep += pyx / self._N return ep def _convergence(self): # 判断是否全部收敛 for last, now in zip(self._lastw, self._w): if abs(last - now) >= self._EPS: return False return True def predict(self, X): # 计算预测概率 Z = self._Zx(X) result = {} for y in self._Y: ss = 0 for x in X: if (x, y) in self._numXY: ss += self._w[self._xyID[(x, y)]] pyx = math.exp(ss) / Z result[y] = pyx return result def train(self, maxiter=1000): # 训练数据 for loop in range(maxiter): # 最大训练次数 print("iter:%d" % loop) self._lastw = self._w[:] for i in range(self._n): ep = self._model_ep(i) # 计算第i个特征的模型期望 self._w[i] += math.log(self._Ep_[i] / ep) / self._C # 更新参数 print("w:", self._w) if self._convergence(): # 判断是否收敛 break
dataset = [['no', 'sunny', 'hot', 'high', 'FALSE'], ['no', 'sunny', 'hot', 'high', 'TRUE'], ['yes', 'overcast', 'hot', 'high', 'FALSE'], ['yes', 'rainy', 'mild', 'high', 'FALSE'], ['yes', 'rainy', 'cool', 'normal', 'FALSE'], ['no', 'rainy', 'cool', 'normal', 'TRUE'], ['yes', 'overcast', 'cool', 'normal', 'TRUE'], ['no', 'sunny', 'mild', 'high', 'FALSE'], ['yes', 'sunny', 'cool', 'normal', 'FALSE'], ['yes', 'rainy', 'mild', 'normal', 'FALSE'], ['yes', 'sunny', 'mild', 'normal', 'TRUE'], ['yes', 'overcast', 'mild', 'high', 'TRUE'], ['yes', 'overcast', 'hot', 'normal', 'FALSE'], ['no', 'rainy', 'mild', 'high', 'TRUE']] maxent = MaxEntropy() x = ['overcast', 'mild', 'high', 'FALSE'] maxent.loadData(dataset) maxent.train() print('predict:', maxent.predict(x))
predict: {'yes': 0.9999971802186581, 'no': 2.819781341881656e-06}
第6章Logistic回归与最大熵模型-习题
习题6.1
确认Logistic分布属于指数分布族。
解答：
第1步：
首先给出指数分布族的定义：
对于随机变量x ，在给定参数η下，其概率分别满足如下形式：
我们称之为指数分布族。
其中：
x ：可以是标量或者向量，可以是离散值也可以是连续值
η ：自然参数
g ( η ) ：归一化系数
h ( x ) ：x的某个函数
**第2步：**证明伯努利分布属于指数分布族
伯努利分布：φ 是y = 1 的概率，即P ( Y = 1 ) = φ
其中
将y 替换成x ，可得
对比可知，伯努利分布属于指数分布族，其中
第3步：
广义线性模型（GLM）必须满足三个假设：
y∣x;θ∼ExponentialFamily(η)，即假设预测变量y 在给定x ，以θ 为参数的条件概率下，属于以η 作为自然参数的指数分布族；
给定x ，求解出以x 为条件的T ( y )的期望E [ T ( y ) ∣ x ]，即算法输出为h ( x ) = E [ T ( y ) ∣ x ]
满足η = θ T x \ ，即自然参数和输入特征向量x 之间线性相关，关系由θ 决定，仅当η是实数时才有意义，若η 是一个向量，则η i = θ i T x
**第4步：**推导伯努利分布的GLM
已知伯努利分布属于指数分布族，对给定的x , η x，求解期望：
可得到Logistic回归算法，故Logistic分布属于指数分布族，得证。
习题6.2
写出Logistic回归模型学习的梯度下降算法。
解答：
对于Logistic模型：
对数似然函数为
似然函数求偏导，可得
梯度函数为：
Logistic回归模型学习的梯度下降算法：
(1) 取初始值x ( 0 ) ∈ R ，置k = 0
(2) 计算f ( x ( k )
(3) 计算梯度g k = g ( x ( k ) 当∥ g k ∥ < ε \|g_k\| 时，停止迭代，令x ∗ = x ( k ) x^* 否则，求λ k使得f ( x ( k ) + λ k g k ) = max λ ⩾ 0 f ( x ( k ) + λ g k )
(4) 置x ( k + 1 ) = x ( k ) + λ k g k x^{(k+1)}=x^{(k)}计算f ( x ( k + 1 ) ) 当∥ f ( x ( k + 1 ) ) − f ( x ( k ) ) ∥ < ε 或 ∥ x ( k + 1 ) − x ( k ) ∥ < ε时，停止迭代，令x ∗ = x ( k + 1 ) x^*
(5) 否则，置k = k + 1 ，转(3)
%matplotlib inline import numpy as np import time import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from pylab import mpl # 图像显示中文 mpl.rcParams['font.sans-serif'] = ['Microsoft YaHei'] class LogisticRegression: def __init__(self, learn_rate=0.1, max_iter=10000, tol=1e-2): self.learn_rate = learn_rate # 学习率 self.max_iter = max_iter # 迭代次数 self.tol = tol # 迭代停止阈值 self.w = None # 权重 def preprocessing(self, X): """将原始X末尾加上一列，该列数值全部为1""" row = X.shape[0] y = np.ones(row).reshape(row, 1) X_prepro = np.hstack((X, y)) return X_prepro def sigmod(self, x): return 1 / (1 + np.exp(-x)) def fit(self, X_train, y_train): X = self.preprocessing(X_train) y = y_train.T # 初始化权重w self.w = np.array([[0] * X.shape[1]], dtype=np.float) k = 0 for loop in range(self.max_iter): # 计算梯度 z = np.dot(X, self.w.T) grad = X * (y - self.sigmod(z)) grad = grad.sum(axis=0) # 利用梯度的绝对值作为迭代中止的条件 if (np.abs(grad) <= self.tol).all(): break else: # 更新权重w 梯度上升——求极大值 self.w += self.learn_rate * grad k += 1 print("迭代次数：{}次".format(k)) print("最终梯度：{}".format(grad)) print("最终权重：{}".format(self.w[0])) def predict(self, x): p = self.sigmod(np.dot(self.preprocessing(x), self.w.T)) print("Y=1的概率被估计为：{:.2%}".format(p[0][0])) # 调用score时，注释掉 p[np.where(p > 0.5)] = 1 p[np.where(p < 0.5)] = 0 return p def score(self, X, y): y_c = self.predict(X) error_rate = np.sum(np.abs(y_c - y.T)) / y_c.shape[0] return 1 - error_rate def draw(self, X, y): # 分离正负实例点 y = y[0] X_po = X[np.where(y == 1)] X_ne = X[np.where(y == 0)] # 绘制数据集散点图 ax = plt.axes(projection='3d') x_1 = X_po[0, :] y_1 = X_po[1, :] z_1 = X_po[2, :] x_2 = X_ne[0, :] y_2 = X_ne[1, :] z_2 = X_ne[2, :] ax.scatter(x_1, y_1, z_1, c="r", label="正实例") ax.scatter(x_2, y_2, z_2, c="b", label="负实例") ax.legend(loc='best') # 绘制p=0.5的区分平面 x = np.linspace(-3, 3, 3) y = np.linspace(-3, 3, 3) x_3, y_3 = np.meshgrid(x, y) a, b, c, d = self.w[0] z_3 = -(a * x_3 + b * y_3 + d) / c ax.plot_surface(x_3, y_3, z_3, alpha=0.5) # 调节透明度 plt.show()
# 训练数据集 X_train = np.array([[3, 3, 3], [4, 3, 2], [2, 1, 2], [1, 1, 1], [-1, 0, 1], [2, -2, 1]]) y_train = np.array([[1, 1, 1, 0, 0, 0]]) # 构建实例，进行训练 clf = LogisticRegression() clf.fit(X_train, y_train) clf.draw(X_train, y_train)
迭代次数：3232次 最终梯度：[ 0.00144779 0.00046133 0.00490279 -0.00999848] 最终权重：[ 2.96908597 1.60115396 5.04477438 -13.43744079]
习题6.3
写出最大熵模型学习的DFP算法。（关于一般的DFP算法参见附录B）
解答：
第1步：
最大熵模型为：
引入拉格朗日乘子，定义拉格朗日函数：
最优化原始问题为：$
对偶问题为：
Ψ(w)称为对偶函数，同时，其解记作
求L ( P , w ) 对P ( y ∣ x ) 的偏导数，并令偏导数等于0，解得：
其中：
则最大熵模型目标函数表示为
第2步：
最大熵模型的DFP算法
输入：目标函数φ ( w )，梯度g ( w ) = ∇ g ( w ) ，精度要求ε
输出：φ ( w ) 的极小值点w ∗
(1)选定初始点w ( 0 ) ，取G 0
为正定对称矩阵，置k = 0
(2)计算g k = g ( w ( k ) ，若 g k < ε \|g_k\| ，则停止计算，得近似解w ∗ = w ( k ) 否则转(3)
(3)置p k = − G k g k (4)一维搜索：求λ k
(5)置w ( k + 1 ) = w ( k ) + λ k p k
(6)计算g k + 1 = g ( w ( k + 1 ) ) ，若∥ g k + 1 ∥ < ε \|g_{k+1}\| ，则停止计算，得近似解w ∗ = w ( k + 1 )否则，按照迭代式算出G k + 1
(7)置k = k + 1 ，转(3)
