天天看點

周志華《機器學習》課後習題(第三章):線性模型

周志華《機器學習》課後習題(第三章):線性模型
對實數集上的函數,可通過求二階導數來判别:若二階導數在區間上非負,則稱為凸函數;若二階導數在區間上恒大于 0,則稱為嚴格凸函數。原書p54)

對于多元函數,其Hessian matrix為半正定即為凸函數。

周志華《機器學習》課後習題(第三章):線性模型
周志華《機器學習》課後習題(第三章):線性模型

3.3 程式設計實作對率回歸,并給出西瓜資料集3.0α上的結果

https://github.com/han1057578619/MachineLearning_Zhouzhihua_ProblemSets/tree/master/ch3--%E7%BA%BF%E6%80%A7%E6%A8%A1%E5%9E%8B/3.3
'''
與原書不同,原書中一個樣本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 折交叉驗證法和留一法所估計出的對率回歸的錯誤率。

https://github.com/han1057578619/MachineLearning_Zhouzhihua_ProblemSets/tree/master/ch3--%E7%BA%BF%E6%80%A7%E6%A8%A1%E5%9E%8B/3.4
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α 上的結果.

https://github.com/han1057578619/MachineLearning_Zhouzhihua_ProblemSets/tree/master/ch3--%E7%BA%BF%E6%80%A7%E6%A8%A1%E5%9E%8B/3.5
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,無法正确分類。即任意一個分類器将樣本分類錯誤,最終結果依然正确,但如果有兩個以上的分類器錯誤,結果就不一定正确了。這是 

周志華《機器學習》課後習題(第三章):線性模型

的由來。

此外,原論文中提到,一個好的糾錯輸出碼應該滿足兩個條件:

  1. 行分離。任意兩個類别之間的codeword距離應該足夠大。
  2. 列分離。任意兩個分類器 
    周志華《機器學習》課後習題(第三章):線性模型
    的輸出應互相獨立,無關聯。這一點可以通過使分類器 
    周志華《機器學習》課後習題(第三章):線性模型
    編碼與其他分類編碼的海明距離足夠大實作,且與其他分類編碼的反碼的海明距離也足夠大(有點繞。)。

第一點其實就是原書提到的,已經解釋過了,說說第二點:

如果兩個分類器的編碼類似或者完全一緻,很多算法(比如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。

繼續閱讀