Hồi quy Logistic (Logistic Regression)¶

Giới thiệu
- So sánh với hồi quy tuyến tính
- Phân loại hồi quy logistic
Hồi quy logistic nhị phân
Hồi quy logistic đa lớp (Multiclass logistic regression)

Giới thiệu ¶

Hồi quy logistic là một thuật toán được sử dụng để phân loại các quan sát theo các danh mục rời rạc. Thay vì có đầu ra là các giá trị liên tục như thuật toán hồi quy tuyến tính, hồi quy logistic đưa đầu ra qua hàm logistic sigmoid để trả về một giá trị biểu thị xác suất mà quan sát đầu vào thể hiện tính chất của hai hay nhiều hơn các danh mục rời rạc.

Ghi chú

Tuy thường được sử dụng trong các bài toán "phân loại", bản thân hồi quy logistic không phải là một thuật toán phân loại mà là một thuật toán hồi quy. Thuật toán này ước lượng xác suất của các kết quả khả dĩ và kết hợp với một số quy tắc (giá trị ngưỡng) để đưa ra quyết định phân loại.
Tên gọi "Logistic" xuất phát từ hàm lôgit được sử dụng trong thuật toán này để tìm ra mối quan hệ giữa bộ dự đoán tuyến tính và trung bình của hàm xác suất.

So sánh với hồi quy tuyến tính ¶

Giả sử ta có một tập dữ liệu bao gồm thời gian học tập và điểm kiểm tra của học sinh. Hồi Quy Tuyến Tính (Linear Regression) và hồi quy logistic sẽ dự đoán theo các cách thức khác nhau:

Hồi quy tuyến tính có thể giúp ta dự đoán điểm kiểm tra của học sinh theo thang điểm 100. Dự đoán của hồi quy tuyến tính là giá trị liên tục (trong một khoảng giá trị).

Hồi quy logistic có thẻ giúp ta dự đoán liệu học sinh đó sẽ đỗ hay trượt bài kiểm tra. Dự đoán của hồi quy logistic là giá trị rời rạc (các giá trị cụ thể hoặc các danh mục). Ta cũng có thể xem được xác suất làm cơ sở cho cách phân loại từng danh mục của mô hình.

Phân loại hồi quy logistic ¶

Nhị phân (Đỗ/Trượt)

Đa lớp (Chó, Mèo, Thỏ, ...)

Thứ tự (Cao, Trung bình, Thấp)

Hồi quy logistic nhị phân ¶

Giả sử ta được cung cấp dữ liệu về kết quả kỳ thi của học sinh và mục tiêu của ta là dự đoán liệu một học sinh sẽ đỗ hay trượt bài kiểm tra dựa trên số giờ ngủ và số giờ học bài của học sinh đó. Ta có hai đặc trưng (số giờ ngủ, số giờ học) và hai danh mục: đỗ (1) và trượt (0).

Ta có thể biểu diễn dữ liệu này dưới dạng biểu đồ phân tán như sau.

scatter plot of exam scores — Biểu đồ phân tán dữ liệu¶

Kích hoạt Sigmoid ¶

Để có thể ánh xạ các giá trị dự đoán với xác suất, ta sử dụng hàm sigmoid. Hàm này ánh xạ bất cứ dữ liệu thực nào thành một giá trị trong khoảng từ \(0\) đến \(1\). Trong học máy, ta sử dụng hàm sigmoid để ánh xạ các dự đoán với xác suất.

Công thức toán học

\[S(z) = \frac{1} {1 + e^{-z}}\]

Trong đó:

\(s(z)\) = đầu ra trong khoảng từ \(0\) đến \(1\) (giá trị xác suất ước lượng)
\(z\) = đầu vào của hàm (giá trị dự đoán của thuật toán, ví dụ như \(mx + b\))
\(e\) = hằng số Euler, và là cơ số của logarit tự nhiên

Đồ thị

graph of sigmoid function — Đồ thị hàm số Sigmoid¶

Code

def sigmoid(z):
  	return 1.0 / (1 + np.exp(-z))

Ranh giới quyết định (Decision boundary)¶

Hàm dự đoán của ta lúc này sẽ trả về giá trị xác suất trong khoảng từ \(0\) đến \(1\). Để có thể ánh xạ xác suất này tới các danh mục rời rạc (đúng/sai, chó/mèo), ta cần chọn một giá trị ngưỡng mà nếu xác suất lớn hơn giá trị này thì ta sẽ phân loại thành danh mục đó, còn nếu thấp hơn thì ta phân loại thành danh mục còn lại.

\[\begin{split}p \geq 0.5, class=1 \\ p < 0.5, class=0\end{split}\]

Ví dụ, nếu giá trị ngưỡng là \(0.5\) và hàm dự đoán trả về \(0.7\), ta có thể phân loại điểm dữ liệu đó là dương tính. Nếu dự đoán là \(0.2\) thì ta có thể phân loại điểm đữ liệu đó là âm tính. Đối với hồi quy tuyến tính đa lớp, ta có thể chọn danh mục mà có xác suất dự đoán cao nhất.

decision boundary — Ranh giới quyết định (màu xanh lam). Các giá trị ở phía trên ranh giới là dương tính, còn lại là âm tính.¶

Đưa ra dự đoán ¶

Với các kiến thức về hàm sigmoid và ranh giới quyết định, lúc này ta có thể viết hàm dự đoán cho mô hình. Một hàm dự đoán trong hồi quy logistic trả về xác suất mà quan sát đầu vào là dương tính, nghĩa là việc quan sát đầu vào thuộc danh mục đó là Đúng (True). Ta gọi danh mục này là \(1\) và được ký hiệu là \(P(class=1)\). Khi xác suất này tiến gần tới \(1\), mô hình sẽ càng chắc chắn hơn là quan sát đó thuộc danh mục \(1\).

Ở phần này, ta sử dụng phương trình dự đoán hồi quy tuyến tính đa biến từ bài hướng dẫn trước về hồi quy tuyến tính.

\[z = W_0 + W_1 \times \text{\{Số giờ học\}} + W_2 \times \text{\{Số giờ ngủ\}}\]

Khác với bài hồi quy tuyến tính, ở đây ta sẽ biến đổi đầu ra qua hàm sigmoid để trả về giá trị xác suất từ \(0\) đến \(1\).

\[P(class=1) = \frac{1} {1 + e^{-z}}\]

Nếu mô hình trả về \(0.4\) thì có nghĩa là mô hình tin rằng học sinh đó chỉ có 40% cơ hội đỗ. Nếu ranh giới quyết định bằng \(0.5\), ta có thể phân loại quan sát này là "Trượt".

Code

Ta sẽ lồng hàm sigmoid ra ngoài hàm dự đoán của hồi quy tuyến tính đa biến.

def predict(features, weights):
    '''
    Trả về một mảng 1D là các xác suất
    mà nhãn danh mục đó == 1
    '''
    z = np.dot(features, weights)
    return sigmoid(z)

Hàm chi phí ¶

Đáng tiếc rằng ta không thể (hay đúng hơn là không nên) sử dụng cùng hàm chi phí trung bình bình phương sai số như trong hồi quy tuyến tính. Chương 3 trong cuốn sách về học sâu (deep learning) của Michael Neilson 5 giải thích rất chi tiết về vấn đề này về mặt toán học, tuy nhiên ở đây ta chỉ cần hiểu đơn giản rằng lý do là vì hàm dự đoán trong hồi quy logistic là phi tuyến (do biến đổi sigmoid). Bình phương dự đoán như trong MSE sẽ tạo ra một hàm không lồi (non-convex) với rất nhiều cực tiểu cục bộ. Nếu hàm chi phí có nhiều cực tiểu cục bộ, thuật toán hạ gradient sẽ khó mà có thể tìm được điểm tối ưu tại cực tiểu toàn cục.

Công thức toán học

Thay vì sử dụng MAE, ta sẽ sử dụng một hàm chi phí có tên gọi là Entropy chéo - Cross-Entropy, hay cũng có thể gọi là hàm mất mát Log. Hàm chi phí entropy chéo được chia làm hai hàm chi phí riêng biệt: một hàm cho \(y=1\) và một cho \(y=0\).

\[\begin{split}J(\theta) = \frac{1}{m} \sum_{i=1}^{m} Cost(h_{\theta}(x^{(i)}), y^{(i)}) \\\end{split}\]

trong đó:

\[\begin{split}& Cost(h_{\theta}(x^{(i)}), y^{(i)}) = -\log(h_{\theta}(x)) && \text{ nếu } y=1 \\ & Cost(h_{\theta}(x^{(i)}), y^{(i)}) = -\log(1 - h_{\theta}(x)) && \text{ nếu } y=0\end{split}\]

Lợi ích của việc sử dụng logarit có thể thấy rõ qua đồ thị của hàm chi với lần lượt với \(y=1\) và \(y=0\). Các hàm số trơn đơn điệu (luôn tăng hoặc luôn giảm) giúp việc tính toán gradient và tối thiểu chi phí dễ dàng hơn.

_images/y1andy2_logistic_function.png — Hình vẽ minh hoạ được trích từ slide về hồi quy logistic của Andrew Ng 1.¶

Ý tưởng chủ đạo của hàm chi phí cần chú ý là nó "phạt" các dự đoán sai với độ chắc chắn cao nhiều hơn so với "thưởng" các dự đoán đúng với độ chắc chắn cao. Cụ thể, ví dụ với \(y=1\), khi mô hình dự đoán đúng và trả về \(h_{\theta}(x)\) tiến tới \(1\) thì hàm entropy chéo trả về chi phí xấp xỉ \(0\), tuy nhiên khi mô hình dự đoán sai với \(h_{\theta}(x)\) tiến tới \(0\) thì giá trị hàm entropy chéo trả về tiến tới vô cực theo cấp số mũ. Hệ quả là với độ chính xác càng cao thì chi phí hàm trả về càng thấp do bản chất logistic của hàm chi phí.

Công thức ngắn gọn

\[J(\theta) = -\frac{1}{m} \sum_{i=1}^{m} [y^{(i)} \log(h_{\theta}(x^{(i)})) + (1-y^{(i)}) \log(1 - h_{\theta}(x^{(i)}))]\]

Phương trình trên nhân từng hàm chi phí với \(y\) và \((1-y)\) để ta có thể sử dụng cùng một hàm cho cả hai trường hợp \(y=1\) và \(y=0\). Nếu \(y=0\), số hạng đầu tiên trong phép cộng bị khử. Nếu \(y=1\), số hạng thứ hai bị khử. Đối với cả hai trường hợp, ta chỉ tính toán đúng vế mà ta cần tính.

Hàm chi phí được vector hoá

\[J(\theta) = \frac{1}{m} \cdot (-y^T log(h) - (1-y)^T log(1-h))\]

trong đó:

\[h = g(X \cdot \theta)\]

Code

def cost_function(features, labels, weights):
    '''
    Sử dụng MAE (Mean Absolute Error)

    Features:(100,3)
    Labels: (100,1)
    Weights:(3,1)

    Trả về mảng 1D các dự đoán
    cost = (labels*log(predictions) + (1-labels)*log(1-predictions) ) / len(labels)
    '''

    # Số quan sát trong tập dữ liệu
    observations = len(labels)

    predictions = predict(features, weights)

    # Tính sai số nếu nhãn = 1
    class1_cost = -labels*np.log(predictions)

    # Tính sai số nếu nhãn = 1
    class2_cost = (1-labels)*np.log(1-predictions)

    # Tính tổng của cả 2 loại sai số
    cost = class1_cost - class2_cost

    # Tính trung bình lỗi làm chi phí
    cost = cost.sum() / observations

    return cost

Hạ gradient ¶

Để tối thiểu hoá mất mát, ta sử dụng Hạ Gradient (Gradient Descent) giống như mô tả trong Hồi Quy Tuyến Tính (Linear Regression). Có nhiều thuật toán tối ưu khác phức tạp hơn ví dụ như gradient liên hợp, tuy nhiên bạn không cần phải quá tập trung vào các thuật toán này. Các thư viện học máy như Scikit-learn hỗ trợ toàn bộ phần lập trình chúng để bạn có thể tập trung vào các vấn đề khác thú vị hơn.

Công thức toán học

Một trong những tính chất hay nhất của hàm sigmoid là công thức đạo hàm của nó rất đẹp. Nếu bạn quan tâm, có một bài khá hay mô tả chi tiết đạo hàm này trên stackoverflow 6. Tác giả Michael Neilson cũng có trình bày vấn đề này trong chương 3 trong quyển sách về học sâu của ông ấy.

\[s'(z) = s(z)(1 - s(z))\]

từ đây ta có thể suy ra công thức đạo hàm rất đẹp và dễ tính toán của hàm chi phí:

\[C' = x(s(z) - y)\]

trong đó:

\(C'\) là đạo hàm chi phí theo các trọng số.

\(y\) là nhãn thực mô tả danh mục quan sát được (\(0\) hoặc \(1\)).

\(s(z)\) là dự đoán của mô hình.

\(x\) là đặc trưng, hay vector các đặc trưng.

Chú ý rằng gradient này giống với gradient của MSE, khác biệt duy nhất nằm ở hàm giả định (hypothesis function).

Mã giả

Repeat {

  1. Tính trung bình gradient
  2. Nhân với tốc độ học
  3. Trừ tích trên khỏi các trọng số

}

Code

def update_weights(features, labels, weights, lr):
    '''
    Thuật toán Hạ Gradient được vector hoá

    Features:(200, 3)
    Labels: (200, 1)
    Weights:(3, 1)
    '''
    N = len(features)

    #1 - Dự đoán kết quả
    predictions = predict(features, weights)

    #2 -  Chuyển vị ma trận đặc trưng từ kích thước (200, 3) về (3, 200)
    # để ta có thể nhân với ma trận sai số (200,1).
    # Trả về một ma trận (3,1) gồm có 3 đạo hàm riêng -
    # mỗi đạo hàm cho một đặc trưng -- đại diện cho tổng
    # độ nghiêng của hàm chi phí qua tất cả các quan sát.
    gradient = np.dot(features.T,  predictions - labels)

    #3 - Tính trung bình đạo hàm của hàm chi phí với mỗi đặc trưng
    gradient /= N

    #4 - Nhân gradient với tốc độ học
    gradient *= lr

    #5 - Cập nhật trọng số bằng cách trừ đi gradient để tối thiểu hoá chi phí
    weights -= gradient

    return weights

Ghép xác suất với các lớp ¶

Bước cuối cùng là gán nhãn các danh mục (\(0\) hoặc \(1\)) cho các giá trị xác suất mô hình dự đoán.

Ranh giới quyết định

def decision_boundary(prob):
  return 1 if prob >= .5 else 0

Ánh xạ xác suất sang các danh mục

def classify(predictions):
  '''
  input  - N element array of predictions between 0 and 1
  output - N element array of 0s (False) and 1s (True)
  '''
  decision_boundary = np.vectorize(decision_boundary)
  return decision_boundary(predictions).flatten()

Ví dụ đầu ra

Probabilities = [ 0.967, 0.448, 0.015, 0.780, 0.978, 0.004]
Classifications = [1, 0, 0, 1, 1, 0]

Huấn luyện ¶

Code huấn luyện mô hình trong ví dụ này giống với đoạn code sử dụng trong hồi quy tuyến tính.

def train(features, labels, weights, lr, iters):
    cost_history = []

    for i in range(iters):
        weights = update_weights(features, labels, weights, lr)

        #Calculate error for auditing purposes
        cost = cost_function(features, labels, weights)
        cost_history.append(cost)

        # Log Progress
        if i % 1000 == 0:
            print("iter: "+str(i) + " cost: "+str(cost))

    return weights, cost_history

Đánh giá mô hình ¶

Nếu mô hình của chúng ta thực sự hoạt động, ta sẽ thấy chi phí giảm dần sau mỗi vòng lặp.

iter: 0 cost: 0.635
iter: 1000 cost: 0.302
iter: 2000 cost: 0.264

Chi phí sau huấn luyện: 0.2487. Trọng số sau huấn luyện: [-8.197, .921, .738]

Nhật ký huấn luyện

_images/logistic_regression_loss_history.png — Đồ thị chi phí của hồi quy logistic trong quá trình huấn luyện.¶

Độ chính xác

Độ chính xác đo khả năng dự đoán của mô hình. Trong trường hợp này, ta chỉ đơn giản so sánh nhãn được dự đoán với nhãn thực tế và chia cho tổng số các quan sát kiểm tra.

def accuracy(predicted_labels, actual_labels):
    diff = predicted_labels - actual_labels
    return 1.0 - (float(np.count_nonzero(diff)) / len(diff))

Ranh giới quyết định

Một kỹ thuật hữu ích khác để có thể biểu diễn độ chính xác của mô hình một cách trực quan là vẽ ranh giới quyết định trên cùng đồ thị với các dự đoán của mô hình để so sánh nhãn dự đoán với nhãn thực. Để thực hiện việc này, ta cần vẽ đồ thị các xác suất mà mô hình dự đoán và tô màu các quan sát dựa theo giá trị nhãn thực. Trong đồ thị dưới, các xác suất phía trên ranh giới quyết định sẽ được phân loại là "True", tuy nhiên mô hình phân loại sai một số quan sát thành "False" (màu đỏ), và ngược lại đối với vùng phía dưới ranh giới.

_images/logistic_regression_final_decision_boundary.png — Ranh giới quyết định sau huấn luyện.¶

Code để vẽ đồ thị với ranh giới quyết định

def plot_decision_boundary(trues, falses):
    fig = plt.figure()
    ax = fig.add_subplot(111)

    no_of_preds = len(trues) + len(falses)

    ax.scatter([i for i in range(len(trues))], trues, s=25, c='b', marker="o", label='Trues')
    ax.scatter([i for i in range(len(falses))], falses, s=25, c='r', marker="s", label='Falses')

    plt.legend(loc='upper right');
    ax.set_title("Decision Boundary")
    ax.set_xlabel('N/2')
    ax.set_ylabel('Predicted Probability')
    plt.axhline(.5, color='black')
    plt.show()

Hồi quy logistic đa lớp (Multiclass logistic regression)¶

Thay vì \(y = {0,1}\), ta sẽ mở rộng định nghĩa này thành \(y = {0,1...n}\). Ý tưởng cơ bản để giải quyết bài toàn này là ta sẽ chạy thuật toán phân loại nhị phân nhiều lần, mỗi lần tương ứng với phân loại một danh mục.

Thuật toán ¶

Chia nhỏ bài toán thành \(n\) bài toán phân loại nhị phân.
Với mỗi danh mục:
- Dự đoán xác suất mà các quan sát thuộc danh mục đó.
\(\text{dự đoán} = max(\text{xác suất quan sát thuộc các danh mục})\)

Với mỗi bài toán con, ta chọn một danh mục (\(1\)) và coi tất cả các danh mục còn lại là một danh mục duy nhất (\(0\)). Sau khi hoàn thành dự đoán xác suất của tất cả các bài toán con, ta phân loại quan sát vào danh mục với giá trị dự đoán cao nhất.

Kích hoạt Softmax ¶

Hàm kích hoạt softmax (có tên gọi khác là softargmax hay hàm trung bình mũ) là một hàm nhận đầu vào là vector gồm K số thực, và chuẩn hoá nó thành một phân phối xác suất gồm có K xác suất được chia tỉ lệ theo hàm mũ tự nhiên của các số đầu vào. Nghĩa là trước khi áp dụng hàm softmax, một số phần tử của vector có thể mang dấu âm, hoặc lớn hơn \(1\), hay tổng các phần tử vector đầu vào có thể không bằng \(1\); tuy nhiên sau khi áp dụng hàm softmax, mỗi phân tử xác suất đầu ra sẽ có giá trị trong khoảng \([0,1]\), và tổng các xác suất sẽ bằng \(1\) để đảm bảo điều kiện tổng các xác suất của một quan sát bằng \(1\).

Hàm softmax chuẩn được định nghĩa bằng công thức sau.

\[\begin{split}\sigma(z_i) & = \frac{e^{z_{(i)}}}{\sum_{j=1}^K e^{z_{(j)}}} \\ \text{for } & i=1,.,.,.,K \\ \text{and } & z=z_1,.,.,.,z_K\end{split}\]

Nói cách khác: ta áp dụng hàm mũ tự nhiên lên từng phần tử \(z_i\) của vector đầu vào \(z\) và chuẩn hoá các giá trị này bằng cách chia cho tổng tất cả các hàm mũ; việc chuẩn hoá giúp đảm bảo rằng tổng của các phần tử của vector xác suất đầu ra \(\sigma(z)\) bằng \(1\) 9.

Ví dụ với Scikit-Learn ¶

Hãy cùng so sánh hiệu năng mô hình trên với mô hình LogisticRegression cung cấp bởi scikit-learn 8.

import sklearn
from sklearn.linear_model import LogisticRegression
from sklearn.cross_validation import train_test_split

# Normalize grades to values between 0 and 1 for more efficient computation
normalized_range = sklearn.preprocessing.MinMaxScaler(feature_range=(-1,1))

# Extract Features + Labels
labels.shape =  (100,) #scikit expects this
features = normalized_range.fit_transform(features)

# Create Test/Train
features_train,features_test,labels_train,labels_test = train_test_split(features,labels,test_size=0.4)

# Scikit Logistic Regression
scikit_log_reg = LogisticRegression()
scikit_log_reg.fit(features_train,labels_train)

#Score is Mean Accuracy
scikit_score = clf.score(features_test,labels_test)
print 'Scikit score: ', scikit_score

#Our Mean Accuracy
observations, features, labels, weights = run()
probabilities = predict(features, weights).flatten()
classifications = classifier(probabilities)
our_acc = accuracy(classifications,labels.flatten())
print 'Our score: ',our_acc

Độ chính xác của Scikit: 0.88. Độ chính xác mô hình trên: 0.89

Tài liệu tham khảo

1: http://www.holehouse.org/mlclass/06_Logistic_Regression.html
2: http://machinelearningmastery.com/logistic-regression-tutorial-for-machine-learning
3: https://scilab.io/machine-learning-logistic-regression-tutorial/
4: https://github.com/perborgen/LogisticRegression/blob/master/logistic.py
5: http://neuralnetworksanddeeplearning.com/chap3.html
6: http://math.stackexchange.com/questions/78575/derivative-of-sigmoid-function-sigma-x-frac11e-x
7: https://en.wikipedia.org/wiki/Monotoniconotonic_function
8: http://scikit-learn.org/stable/modules/linear_model.html#logistic-regression>
9: https://en.wikipedia.org/wiki/Softmax_function