对实数集上的函数，可通过求二阶导数来判别：若二阶导数在区间上非负，则称为凸函数；若二阶导数在区间上恒大于 0，则称为严格凸函数。原书p54）
对于多元函数，其Hessian matrix为半正定即为凸函数。
3.3 编程实现对率回归，并给出西瓜数据集3.0α上的结果
''' 与原书不同，原书中一个样本xi 为列向量，本代码中一个样本xi为行向量 尝试了两种优化方法，梯度下降和牛顿法。两者结果基本相同，不过有时因初始化的原因， 会导致牛顿法中海森矩阵为奇异矩阵，np.linalg.inv(hess)会报错。以后有机会再写拟牛顿法吧。 ''' import numpy as np import pandas as pd from matplotlib import pyplot as plt from sklearn import linear_model def sigmoid(x): s = 1 / (1 + np.exp(-x)) return s def J_cost(X, y, beta): ''' :param X: sample array, shape(n_samples, n_features) :param y: array-like, shape (n_samples,) :param beta: the beta in formula 3.27 , shape(n_features + 1, ) or (n_features + 1, 1) :return: the result of formula 3.27 ''' X_hat = np.c_[X, np.ones((X.shape[0], 1))] beta = beta.reshape(-1, 1) y = y.reshape(-1, 1) Lbeta = -y * np.dot(X_hat, beta) + np.log(1 + np.exp(np.dot(X_hat, beta))) return Lbeta.sum() def gradient(X, y, beta): ''' compute the first derivative of J(i.e. formula 3.27) with respect to beta i.e. formula 3.30 ---------------------------------- :param X: sample array, shape(n_samples, n_features) :param y: array-like, shape (n_samples,) :param beta: the beta in formula 3.27 , shape(n_features + 1, ) or (n_features + 1, 1) :return: ''' X_hat = np.c_[X, np.ones((X.shape[0], 1))] beta = beta.reshape(-1, 1) y = y.reshape(-1, 1) p1 = sigmoid(np.dot(X_hat, beta)) gra = (-X_hat * (y - p1)).sum(0) return gra.reshape(-1, 1) def hessian(X, y, beta): ''' compute the second derivative of J(i.e. formula 3.27) with respect to beta i.e. formula 3.31 ---------------------------------- :param X: sample array, shape(n_samples, n_features) :param y: array-like, shape (n_samples,) :param beta: the beta in formula 3.27 , shape(n_features + 1, ) or (n_features + 1, 1) :return: ''' X_hat = np.c_[X, np.ones((X.shape[0], 1))] beta = beta.reshape(-1, 1) y = y.reshape(-1, 1) p1 = sigmoid(np.dot(X_hat, beta)) m, n = X.shape P = np.eye(m) * p1 * (1 - p1) assert P.shape[0] == P.shape[1] return np.dot(np.dot(X_hat.T, P), X_hat) def update_parameters_gradDesc(X, y, beta, learning_rate, num_iterations, print_cost): ''' update parameters with gradient descent method -------------------------------------------- :param beta: :param grad: :param learning_rate: :return: ''' for i in range(num_iterations): grad = gradient(X, y, beta) beta = beta - learning_rate * grad if (i % 10 == 0) & print_cost: print('{}th iteration, cost is {}'.format(i, J_cost(X, y, beta))) return beta def update_parameters_newton(X, y, beta, num_iterations, print_cost): ''' update parameters with Newton method :param beta: :param grad: :param hess: :return: ''' for i in range(num_iterations): grad = gradient(X, y, beta) hess = hessian(X, y, beta) beta = beta - np.dot(np.linalg.inv(hess), grad) if (i % 10 == 0) & print_cost: print('{}th iteration, cost is {}'.format(i, J_cost(X, y, beta))) return beta def initialize_beta(n): beta = np.random.randn(n + 1, 1) * 0.5 + 1 return beta def logistic_model(X, y, num_iterations=100, learning_rate=1.2, print_cost=False, method='gradDesc'): ''' :param X: :param y:~ :param num_iterations: :param learning_rate: :param print_cost: :param method: str 'gradDesc' or 'Newton' :return: ''' m, n = X.shape beta = initialize_beta(n) if method == 'gradDesc': return update_parameters_gradDesc(X, y, beta, learning_rate, num_iterations, print_cost) elif method == 'Newton': return update_parameters_newton(X, y, beta, num_iterations, print_cost) else: raise ValueError('Unknown solver %s' % method) def predict(X, beta): X_hat = np.c_[X, np.ones((X.shape[0], 1))] p1 = sigmoid(np.dot(X_hat, beta)) p1[p1 >= 0.5] = 1 p1[p1 < 0.5] = 0 return p1 if __name__ == '__main__': data_path = r'C:\Users\hanmi\Documents\xiguabook\watermelon3_0_Ch.csv' # data = pd.read_csv(data_path).values is_good = data[:, 9] == '是' is_bad = data[:, 9] == '否' X = data[:, 7:9].astype(float) y = data[:, 9] y[y == '是'] = 1 y[y == '否'] = 0 y = y.astype(int) plt.scatter(data[:, 7][is_good], data[:, 8][is_good], c='k', marker='o') plt.scatter(data[:, 7][is_bad], data[:, 8][is_bad], c='r', marker='x') plt.xlabel('密度') plt.ylabel('含糖量') # 可视化模型结果 beta = logistic_model(X, y, print_cost=True, method='gradDesc', learning_rate=0.3, num_iterations=1000) w1, w2, intercept = beta x1 = np.linspace(0, 1) y1 = -(w1 * x1 + intercept) / w2 ax1, = plt.plot(x1, y1, label=r'my_logistic_gradDesc') lr = linear_model.LogisticRegression(solver='lbfgs', C=1000) # 注意sklearn的逻辑回归中，C越大表示正则化程度越低。 lr.fit(X, y) lr_beta = np.c_[lr.coef_, lr.intercept_] print(J_cost(X, y, lr_beta)) # 可视化sklearn LogisticRegression 模型结果 w1_sk, w2_sk = lr.coef_[0, :] x2 = np.linspace(0, 1) y2 = -(w1_sk * x2 + lr.intercept_) / w2 ax2, = plt.plot(x2, y2, label=r'sklearn_logistic') plt.legend(loc='upper right') plt.show()
3.4 选择两个 UCI 数据集，比较 10 折交叉验证法和留一法所估计出的对率回归的错误率。
import numpy as np from sklearn import linear_model from sklearn.model_selection import LeaveOneOut from sklearn.model_selection import cross_val_score data_path = r'C:\Users\hanmi\Documents\xiguabook\Transfusion.txt' data = np.loadtxt(data_path, delimiter=',').astype(int) X = data[:, :4] y = data[:, 4] m, n = X.shape # normalization X = (X - X.mean(0)) / X.std(0) # shuffle index = np.arange(m) np.random.shuffle(index) X = X[index] y = y[index] # 使用sklarn 中自带的api先 # k-10 cross validation lr = linear_model.LogisticRegression(C=2) score = cross_val_score(lr, X, y, cv=10) print(score.mean()) # LOO loo = LeaveOneOut() accuracy = 0 for train, test in loo.split(X, y): lr_ = linear_model.LogisticRegression(C=2) X_train = X[train] X_test = X[test] y_train = y[train] y_test = y[test] lr_.fit(X_train, y_train) accuracy += lr_.score(X_test, y_test) print(accuracy / m) # 两者结果几乎一样 # 自己写一个试试 # k-10 # 这里就没考虑最后几个样本了。 num_split = int(m / 10) score_my = [] for i in range(10): lr_ = linear_model.LogisticRegression(C=2) test_index = range(i * num_split, (i + 1) * num_split) X_test_ = X[test_index] y_test_ = y[test_index] X_train_ = np.delete(X, test_index, axis=0) y_train_ = np.delete(y, test_index, axis=0) lr_.fit(X_train_, y_train_) score_my.append(lr_.score(X_test_, y_test_)) print(np.mean(score_my)) # LOO score_my_loo = [] for i in range(m): lr_ = linear_model.LogisticRegression(C=2) X_test_ = X[i, :] y_test_ = y[i] X_train_ = np.delete(X, i, axis=0) y_train_ = np.delete(y, i, axis=0) lr_.fit(X_train_, y_train_) score_my_loo.append(int(lr_.predict(X_test_.reshape(1, -1)) == y_test_)) print(np.mean(score_my_loo)) # 结果都是类似
3.5 编辑实现线性判别分析，并给出西瓜数据集 3.0α 上的结果.
import numpy as np import pandas as pd from matplotlib import pyplot as plt class LDA(object): def fit(self, X_, y_, plot_=False): pos = y_ == 1 neg = y_ == 0 X0 = X_[neg] X1 = X_[pos] u0 = X0.mean(0, keepdims=True) # (1, n) u1 = X1.mean(0, keepdims=True) sw = np.dot((X0 - u0).T, X0 - u0) + np.dot((X1 - u1).T, X1 - u1) w = np.dot(np.linalg.inv(sw), (u0 - u1).T).reshape(1, -1) # (1, n) if plot_: fig, ax = plt.subplots() ax.spines['right'].set_color('none') ax.spines['top'].set_color('none') ax.spines['left'].set_position(('data', 0)) ax.spines['bottom'].set_position(('data', 0)) plt.scatter(X1[:, 0], X1[:, 1], c='k', marker='o', label='good') plt.scatter(X0[:, 0], X0[:, 1], c='r', marker='x', label='bad') plt.xlabel('密度', labelpad=1) plt.ylabel('含糖量') plt.legend(loc='upper right') x_tmp = np.linspace(-0.05, 0.15) y_tmp = x_tmp * w[0, 1] / w[0, 0] plt.plot(x_tmp, y_tmp, '#808080', linewidth=1) wu = w / np.linalg.norm(w) # 正负样板店 X0_project = np.dot(X0, np.dot(wu.T, wu)) plt.scatter(X0_project[:, 0], X0_project[:, 1], c='r', s=15) for i in range(X0.shape[0]): plt.plot([X0[i, 0], X0_project[i, 0]], [X0[i, 1], X0_project[i, 1]], '--r', linewidth=1) X1_project = np.dot(X1, np.dot(wu.T, wu)) plt.scatter(X1_project[:, 0], X1_project[:, 1], c='k', s=15) for i in range(X1.shape[0]): plt.plot([X1[i, 0], X1_project[i, 0]], [X1[i, 1], X1_project[i, 1]], '--k', linewidth=1) # 中心点的投影 u0_project = np.dot(u0, np.dot(wu.T, wu)) plt.scatter(u0_project[:, 0], u0_project[:, 1], c='#FF4500', s=60) u1_project = np.dot(u1, np.dot(wu.T, wu)) plt.scatter(u1_project[:, 0], u1_project[:, 1], c='#696969', s=60) ax.annotate(r'u0 投影点', xy=(u0_project[:, 0], u0_project[:, 1]), xytext=(u0_project[:, 0] - 0.2, u0_project[:, 1] - 0.1), size=13, va="center", ha="left", arrowprops=dict(arrowstyle="->", color="k", ) ) ax.annotate(r'u1 投影点', xy=(u1_project[:, 0], u1_project[:, 1]), xytext=(u1_project[:, 0] - 0.1, u1_project[:, 1] + 0.1), size=13, va="center", ha="left", arrowprops=dict(arrowstyle="->", color="k", ) ) plt.axis("equal") # 两坐标轴的单位刻度长度保存一致 plt.show() self.w = w self.u0 = u0 self.u1 = u1 return self def predict(self, X): project = np.dot(X, self.w.T) wu0 = np.dot(self.w, self.u0.T) wu1 = np.dot(self.w, self.u1.T) return (np.abs(project - wu1) < np.abs(project - wu0)).astype(int) if __name__ == '__main__': data_path = r'C:\Users\hanmi\Documents\xiguabook\watermelon3_0_Ch.csv' data = pd.read_csv(data_path).values X = data[:, 7:9].astype(float) y = data[:, 9] y[y == '是'] = 1 y[y == '否'] = 0 y = y.astype(int) lda = LDA() lda.fit(X, y, plot_=True) print(lda.predict(X)) # 和逻辑回归的结果一致 print(y)
3.6 线性判别分析仅在线性可分数据上能获得理想结果，试设计一个改进方法，使其能较好地周于非线性可分数据。
答：
引入核函数，原书p137，有关于核线性判别分析的介绍。
3.7 令码长为 9，类别数为 4，试给出海明距离意义下理论最优的 ECOC二元码并证明之。
答：
原书对很多地方解释没有解释清楚，把原论文看了一下《Solving Multiclass Learning Problems via Error-Correcting Output Codes》。
先把几个涉及到的理论解释一下。
首先原书中提到：
对同等长度的编码，理论上来说，任意两个类别之间的编码距离越远，则纠错能力越强。因此，在码长较小时可根据这个原则计算出理论最优编码。
拿上图论文中的例子解释一下，上图中，所有类别之间的海明距离都为4，假设一个样本正确的类别为 
，那么codeword应该为 ‘0 0 1 1 0 0 1 1’，若此时有一个分类器输出错误，变成‘0 0 0 1 0 0 1 1’，那么此时距离最近的仍然为 
，若有两个分类输出错误如‘0 0 0 0 0 0 1 1’，此时与 
的海明距离都为2，无法正确分类。即任意一个分类器将样本分类错误，最终结果依然正确，但如果有两个以上的分类器错误，结果就不一定正确了。这是 
的由来。
此外，原论文中提到，一个好的纠错输出码应该满足两个条件：
- 行分离。任意两个类别之间的codeword距离应该足够大。
- 列分离。任意两个分类器
的输出应相互独立，无关联。这一点可以通过使分类器 
编码与其他分类编码的海明距离足够大实现，且与其他分类编码的反码的海明距离也足够大（有点绕。）。
第一点其实就是原书提到的，已经解释过了，说说第二点：
如果两个分类器的编码类似或者完全一致，很多算法（比如C4.5）会有相同或者类似的错误分类，如果这种同时发生的错误过多，会导致纠错输出码失效。（翻译原论文）
个人理解就是：若增加两个类似的编码，那么当误分类时，就从原来的1变成3，导致与真实类别的codeword海明距离增长。极端情况，假设增加两个相同的编码，此时任意两个类别之间最小的海明距离不会变化依然为 
，而纠错输出码输出的codeword与真实类别的codeword的海明距离激增（从1变成3）。所以如果有过多同时发出的错误分类，会导致纠错输出码失效。
另外，两个分类器的编码也不应该互为反码，因为很多算法（比如C4.5，逻辑回归）对待0-1分类其实是对称的，即将0-1类互换，最终训练出的模型是一样的。也就是说两个编码互为补码的分类器是会同时犯错的。同样也会导致纠错输出码失效。
当然当类别较少时，很难满足上面这些条件。如上图中，一共有三类，那么只有 
中可能的分类器编码（ 
），其中后四种（ 
）是前四种的反码，都应去除，再去掉全为0的
，就只剩下三种编码选择了，所以很难满足上述的条件。事实上，对于 
种类别的分类，再去除反码和全是0或者1的编码后，就剩下 
中可行的编码。
原论文中给出了构造编码的几种方法。其中一个是：
当码长为9时，那么 
之后加任意两个编码，即为最优编码，因为此时再加任意的编码都是先有编码的反码，此时，类别之间最小的海明距离都为4，不会再增加。
3.8 ECOC 编码能起到理想纠错作用的重要条件是：在每一位编码上出错的概率相当且独立。试析多分类任务经 ECOC 编码后产生的二类分类器满足该条件的可能性及由此产生的影响。
答：
条件分解为两个：一是出错的概率相当，二是出错的可能性相互独立。
先看第一个把，其实就是每个一位上的分类器的泛化误差相同，要满足这个条件其实取决于样本之间的区分难度，若两个类别本身就十分相似，即越难区分，训练出的分类器出错的概率越大，原书p66也提到：
将多个类拆解为两个"类别子集“，所形成的两个类别子集的区分难度往往不同，即其导致的二分类问题的难度不同。
所以每个编码拆解后类别之间的差异越相同（区分难度相当），则满足此条件的可能性越大。在实际中其实很难满足。
第二个，相互独立。在3.7中也提到过，原论文中也提出一个好的纠错输出码应该满足的其中一个条件就是各个位上分类器相互独立，当类别越多时，满足这个条件的可能性越大，在3.7中也解释了当类别较少时，很难满足这个条件。
至于产生的影响。西瓜书上也提到：
一个理论纠错牲质很好、但导致的三分类问题较难的编码，与另一个理论纠错性质差一些、但导致的二分类问题较简单的编码，最终产生的模型性能孰强孰弱很难说。
3.9 使用 OvR 和 MvM 将多分类任务分解为二分类任务求解时，试述为何无需专门针对类别不平衡性进行处理。
答：
p66 其实已经给出答案了：
对 OvR 、 MvM 来说，由于对每个类进行了相同的处理，其拆解出的二分类任务中类别不平衡的影响会相互抵消，因此通常不需专门处理.
3.10 试推导出多分类代价敏感学习(仅考虑基于类别的误分类代价)使用"再缩放"能获得理论最优解的条件。
答：
这道题目其实是周志华教授的一篇论文《On Multi-Class Cost-Sensitive Learning》。把论文理论部分读了一遍。现在尝试概述一遍吧。
首先说一点关于“再缩放”的个人理解：无论是代价敏感学习还是非代价敏感学习中，“再缩放”各种方法（过采样、欠采样、阈值移动等）都是在调整各类别对模型的影响程度，即各类别的权重。
在《On Multi-Class Cost-Sensitive Learning》中，引用了另外一篇论文《The Foundations of Cost-Sensitive Learning》的一个理论：
其伴随矩阵秩小于c。
