1 Bài toán ước lượng xác xuất

Trong thực tế chúng ta thường gặp những bài toán về ước lượng xác xuất để một sự kiện có khả năng xảy ra là bao nhiêu. Đối với những bài toán này thì miền giá trị của biến target chỉ nằm trong khoảng từ [0,1] nên việc áp dụng các model linear regression sẽ không hợp lý do hàm tuyến tính không bị giới hạn miền giá trị trả về. Nếu áp dụng các bài toán dạng non-parametric như K nearest neighbor chúng ta sẽ không ước lượng được xác xuất để tại một data point khả năng để biến cố xảy ra là bao nhiêu. Do đó chúng ta cần xây dựng một hàm xác xuất mới đảm bảo được các tiêu chí:

  • Hàm số phải liên tục và nhận miền giá trị trong khoảng (0,1).

  • Đạo hàm của nó phải là một hàm số đẹp. Khi đó việc áp dụng thuật toán gradient descent sẽ đơn giản hơn.

  • Hàm số phải đồng biến hoặc nghịch biến đối với các biến predictor bởi trong thực tế một biến sẽ chỉ đóng 1 vai trò ảnh hưởng tiêu cực hoặc tích cực lên biến target mà không thể đóng 2 vai trò cùng lúc (vừa tác động cùng chiều và ngược chiều).

  • Hàm loss function phải tìm được nghiệm tối ưu để xác xuất tiên nghiệm đối với một bộ dữ liệu biết trước là lớn nhất. Hay nói cách khác tồn tại ước lượng hợp lý tối đa (maximum loglikelihood estimation) cho hàm loss function.

Dựa trên các giả định trên, hàm logistic đã được tìm ra để đáp ứng được các tiêu chí đó. Bên dưới chúng ta sẽ đi vào một ví dụ cụ thể.

Ví dụ về logistic:

Chúng ta sẽ lấy ví dụ về một mẫu số năm kinh nghiệm làm việc và khả năng làm xếp:

ExpYears <- c(1,2,2.5,6.6,2.8,4,5,1.6,4,2.8,7,8,1.5,3.5,4.5,3.7,9,1.2,5.5,6)
Boss <-     c(0,0,0  ,1  ,1  ,1,1,0  ,1,0  ,0,1,0  ,0  ,0  ,1  ,1,0  ,1  ,1)
(df <- data.frame(ExpYears, Boss))

Trên thực tế chúng ta thấy rằng số năm kinh nghiệm thường sẽ tác động dương đến khả năng làm xếp. Một người có số năm kinh nghiệm làm việc nhiều sẽ có khả năng làm xếp cao hơn so với những người có số năm kinh nghiệm làm việc ít. Tất nhiên cũng có những trường hợp ngoại lệ như làm việc 1-2 năm đã được cất nhắc làm xếp hoặc làm việc 7 năm nhưng vẫn chỉ là nhân viên bình thường. Điều này là do số năm kinh nghiệm không phải là thước đo duy nhất bên cạnh các yếu tố khác như môi trường làm việc, ngoại ngữ, mối quan hệ, khả năng gây ảnh hưởng,….

Do chúng ta đã thấy rằng không phải 100% những nhân viên làm việc > 5 năm sẽ có khả năng làm xếp và không phải nhân viên làm việc từ 1-2 năm sẽ không thể làm xếp nên trong biểu diễn số năm kinh nghiệm thì 2 nhóm 1 và 0 (làm xếp và không làm xếp) sẽ bị overlap lên nhau. Nếu áp dụng thuật toán Perceptron (PLA) sẽ không hợp lý do không tồn tại một đường boundary để phân chia chính xác hoàn toàn các điểm này.

Chúng ta biểu diễn trên đồ thị theo 2 nhóm để thấy rõ hơn:

library(ggplot2)
library(dplyr)
df %>% mutate(Boss = Boss %>% as.character()) %>% ggplot(aes(x = ExpYears, y = Boss)) +
  geom_point(aes(group = Boss, color = Boss, shape = Boss), size = 2) +
  scale_x_continuous(breaks = seq(1,10,1))

Trong thuật toán perceptron sẽ không có trục Boss do nó chỉ biểu diễn các điểm dựa trên tọa độ của biến predictor. Tuy nhiên mình vẽ thêm trục Boss để chúng ta hình dung dẽ hơn 2 nhóm làm xếp và không làm xếp. Trục hoành biểu diễn năm kinh nghiệm. Khi chập 2 đường thẳng y = 1 và y = 0 thì 2 nhóm này overlap nhau ở những điểm mà thuộc cùng range năm kinh nghiệm của cả nhóm 0 và 1.

2 Xây dựng sigmoid function

Giả sử \(\bar{\mathbf{X}} = [\mathbf{x_1},\mathbf{x_2},...,\mathbf{x_N}] \in \mathbb{R}^{(d+1) \times N}\) là matrix mở rộng của \(\mathbf{X}\) chứa các cột là các điểm dữ liệu \(\mathbf{x_i}\) thuộc không gian \(\mathbb{R}^{(d+1) \times 1}\) (\(\mathbf{x_i}\) có phần tử đầu tiên là \(x_{i1} = 1\)). Vector hệ số của phương trình hồi qui \(\mathbf{y}\) theo \(\mathbf{X}\) là \(\mathbf{w}^T = [w_0, w_1, ..., w_d] \in \mathbb{R}^{1 \times (d+1)}\) với \(w_0\) là hệ số tự do. Xác xuất để rơi vào lần lượt các class 0 và 1 tại một điểm dữ liệu \(\mathbf{x_i}\) và vector hệ số \(\mathbf{w}\):

\[P(y_i = 1|\mathbf{x_i}; \mathbf{w}) = f(\mathbf{w}^T\mathbf{x_i})\] \[P(y_i = 0|\mathbf{x_i}; \mathbf{w}) = 1 - f(\mathbf{w}^T\mathbf{x_i})\]

Kí hiệu \(P(y_i = 0|\mathbf{x_i}; \mathbf{w})\) là xác xuất có điều kiện tại điểm \(\mathbf{x_i}\) và \(\mathbf{w}\). Chúng ta gọi \(f(\mathbf{w}^T\mathbf{x_i})\) là hàm mật độ xác xuất (density function) tại điểm \(\mathbf{x_i}\). Đặt \(z_i = f(\mathbf{w}^T\mathbf{x_i})\). Kết hợp cả 2 xác xuất đối với class 0 và 1 trên ta được:

\[\begin{eqnarray}P(y_i|\mathbf{x_i}; \mathbf{w}) &=& f(\mathbf{w}^T\mathbf{x_i})^{y_i}(1 - f(\mathbf{w}^T\mathbf{x_i}))^{(1-y_i)} \\&=& z_i^{y_i}(1-z_i)^{(1-y_i)} \end{eqnarray}\]

Ta sẽ tìm một nghiệm \(\mathbf{w}\) sao cho hàm hợp lý tối đa của hàm số đạt giá trị lớn nhất. Quá trình này gọi là maximum likelihood function. Bài toán tìm nghiệm của hàm hợp lý tối đa rất phổ biến trong các dạng toán hồi qui mà kết quả trả về là xác xuất. Bởi hàm hợp lý tối đa cho ta biết với một họ \(\Theta\) các tham số cho trước thì nghiệm \(\theta\) là bao nhiêu để xác xuất tiên nghiệm đối với một tập dữ liệu nào đó là lớn nhất. Nghiệm của bài toán còn gọi là ước lượng hợp lý tối đa và đối với bài toán này chính là giá trị \(\mathbf{w}\) sao cho xác xuất \(P(\mathbf{y}| \mathbf{X};\mathbf{w})\) đạt giá trị lớn nhất được biểu diễn như sau:

\[\mathbf{w} = \arg \max_{_\mathbf{w}} P(\mathbf{y}| \mathbf{X};\mathbf{w})\] Để cho đơn giản chúng ta cho các điểm dữ liệu là độc lập với nhau. Trên thực tế vẫn tồn tại mối tương quan chẳng hạn như một nhóm những nhân viên làm việc cùng nhau nếu một người trong số họ được thăng chức họ có thể kéo theo những đồng nghiệp khác cùng thăng chức. Khi đó xác xuất đối với toàn bộ các điểm dữ liệu là:

\[\begin{eqnarray} P(\mathbf{y}|\mathbf{X}; \mathbf{w}) &=& \prod_{i=1}^N P(y_i| \mathbf{x}_i; \mathbf{w}) \\ &=& \prod_{i=1}^N z_i^{y_i}(1 - z_i)^{1- y_i} \end{eqnarray}\]

Để thuận tiện cho tính đạo hàm ta sẽ lấy logarit hai vế. Hàm thu được gọi là loglikelihood (log của hàm likelihood) và chúng ta cần tìm nghiệm để hàm này đạt giá trị lớn nhất. Đảo dấu 2 vế để chuyển từ giá trị lớn nhất sang giá trị nhỏ nhất ta sẽ thu được hàm loss function của bài toán như sau:

\[\begin{eqnarray} J(\mathbf{w}) &=& -\log P(\mathbf{y}|\mathbf{X}; \mathbf{w})&=&-\sum_{i=1}^N [y_i \log(z_i) + (1- y_i)\log(1 - z_i)] \end{eqnarray}\]

Đạo hàm của hàm loss fucntion tại điểm dữ liệu \(\mathbf{x_i}\) là: \[\begin{eqnarray} \frac{\partial J(\mathbf{w}|\mathbf{x_i}; y_i)}{\partial \mathbf{w}} &=& -(\frac{y_i} {z_i} - \frac{(1- y_i)}{(1 - z_i)})\frac{\partial z_i}{\partial{\mathbf{w}}} \\ &=& \frac{z_i - y_i}{zi(1-z_i)} \frac{\partial z_i}{\partial{\mathbf{w}}} \space (1) \end{eqnarray}\]

Mục đích của ta là tìm một hàm \(z\) sao cho khi tính đạo hàm theo \(\mathbf{w}\) thì mẫu được triệt tiêu để hàm số trở nên gọn hơn. Đặt \(s = \mathbf{w}^T\mathbf{x}\) ta suy ra:

\[ \frac{\partial z_i}{\partial{\mathbf{w}}} = \frac{\partial z_i}{\partial{s}} \frac{\partial s}{\partial{\mathbf{w}}} = \frac{\partial z_i}{\partial{s}} \mathbf{x} \space (2)\]

Để triệt tiêu mẫu thì ta phải chọn \(s\) sao cho: \[\begin{eqnarray}\frac{\partial z}{\partial{s}} &=& z(1-z) \leftrightarrow \frac{\delta z}{z(1-z)} &=&\delta s \\ &\leftrightarrow& (\frac{1}{z}+\frac{1}{1-z}) \delta z &=& \delta s \\ &\leftrightarrow& \log {z} - \log {(1-z)} &=& s \\ &\leftrightarrow& \log \frac{z}{(1-z)} &=& s \\ &\leftrightarrow& \frac{z}{(1-z)} &=& e^s \\ &\leftrightarrow& z &=& \frac{e^s}{1+e^s} = \sigma(s)\space (3) \end{eqnarray}\]

Đây chính là phương trình của hàm sigmoid và chắc chúng ta đã hiểu cách mà hàm sigmoid đã được tìm ra như thế nào. Một hàm số khác được sử dụng khá phổ biến có dạng gần như hàm sigmoid là hàm tanh:

\[tanh(s) = \frac{e^{s}-e^{-s}}{e^{s}+e^{-s}}\]

Dễ dàng suy ra mối liên hệ giữa tanh và sigmoid:

\[tanh(s) = sigmoid(e^{2s}) - sigmoid(e^{-2s})\]

Gradient của hàm sigmoid: thế (2) và (3) vào phương trình số (1) ta được:

\[\begin{eqnarray} \frac{\partial J(\mathbf{w}|\mathbf{x_i}; y_i)}{\partial \mathbf{w}} &=& (z_i - y_i)\mathbf{x_i} \end{eqnarray}\]

Do đó hàm số để cập nhật điểm \(\mathbf{w}\) tại \(\mathbf{x_i}\) là: \[\mathbf{w_{t+1}} = \mathbf{w_{t}} + \eta(y_i-z_i)\mathbf{x_i}\]

3 Hàm entrophy và lý thuyết thông tin

Giả sử với một bộ dữ liệu cho trước khi tuân theo một phân phối thực nghiệm, tại một điểm dữ liệu chúng ta thu được xác xuất ước lượng đối với class i là \(q_i\). Thực tế cho thấy xác xuất để rơi vào class i là \(p_i\). Để xác xuất thu được từ phân phối thực nghiệm sát với xác xuất thực tế ta phải xây dựng một hàm số để nghiệm của nó trả về là phân phối sát nhất với xác xuất thực tế. Hàm số này chính là hàm số entrophy có dạng như bên dưới:

\[H(\mathbf{p},\mathbf{q}) = \mathbf{E_p}[-\log \mathbf{q}]\]

Trong trường hợp \(\mathbf{p}, \mathbf{q}\) là rời rạc ta sẽ biểu diễn entrophy như sau:

\[H(\mathbf{p},\mathbf{q}) = -\sum_{i=1}^{C} p_i\log(q_i)\]

Giá trị của entrophy là kì vọng bình quân của giá trị âm logarit xác xuất thực nghiệm có trọng số là xác xuất thực tế. Hàm số này sẽ đánh trọng số cao đối với những class có khả năng xảy ra cao và trọng số thấp hơn với những class có khả năng xảy ra thấp để lượng thông tin thể hiện là lớn nhất. Người ta đã chứng minh được rằng gia trị nhỏ nhất của entrophy đạt được khi hai phân phối thực nghiệm và thực tế có phân bố như nhau tức là \(p_i = q_i\). Đây là bài toán tối ưu lồi có thể giải được bằng phương pháp nhân tử lagrangian.

bài toán: cho một phân phối xác xuất thực tế \(\mathbf{p} = [p_{1}, p_{2}, ..., p_{C}]^{T} \in\mathbb {R}^{C}\) (\(\mathbb{R}^{C}\) là kí hiệu không gian vector số thực C chiều) thỏa mãn \(p_{i} \geq 0\) và một phân phối thực nghiệm là \(\mathbf{q} = [q_{1}, q_{2}, ..., q_{C}]^{T} \in \mathbb{R}^{C}\) với giả định \(q_{i} > 0, \forall i\) và \(\sum_{i=1}^{N} q_{i} = 1\) . Tìm phân phối thực nghiệm \(\mathbf{q}\) để hàm tối ưu sau đây đạt giá trị nhỏ nhất: \[H(\mathbf{p},\mathbf{q}) = -\sum_{i=1}^{C} p_i\log(q_i)\]

Xây dựng hàm nhân tử lagrangian của bài toán dựa trên điều kiện ràng buộc \(\sum_{i=1}^{C}q_i = 1\):

\[\mathcal{L}(q_1, q_2, \dots , q_n, \lambda) = -\sum_{i=1}^{C} p_i\log(q_i) + \lambda(\sum_{i=1}^{C} q_i - 1)\] Hàm số đạt giá trị cực tiểu khi:

\[\Delta_{q_1, q_2, \dots, q_n,\lambda}\mathcal{L}(q_1, q_2, \dots , q_n, \lambda) = 0 \Leftrightarrow \left\{ \begin{matrix} \frac{-p_i}{q_i} + \lambda &=& 0 ~~, \forall i = 1,\dots, n \\ \sum_{i=1}^{C} q_i &=& 1 \end{matrix} \right.\]

Từ phương trình thứ nhất của hệ phương trình ta suy ra \(\frac{p_i}{q_i} = \lambda, \forall i\) .Mặt khác do \(1 = \sum_{i=1}^{C} q_i = \sum_{i=1}^{C} p_i = \sum_{i=1}^{C} \lambda q_i = \lambda \Leftrightarrow \lambda = 1 \Rightarrow p_i = q_i, \forall i\)

Rõ ràng trong bài toán linear regression chúng ta sử dụng norm bậc 2 để tính khoảng cách giữa 2 giá trị ước lượng và thực tế. Đối với bài toán hồi qui dạng xác xuất chúng ta cũng có thể dựa trên norm bậc 2 để tối ưu. Tuy nhiên ta có thể chứng minh rằng sử dụng entrophy và norm bậc 2 nghiệm sẽ cùng hội tụ tại một điểm. Thật vậy, với cùng giả thuyết như trên, hàm loss function trong trường hợp norm bậc 2:

\[J(\mathbf{p},\mathbf{q}) = -\sum_{i=1}^{C} (p_i-q_i)^2\]

Hàm lagrangian:

\[\mathcal{L}(q_1, q_2, \dots , q_n, \lambda) = \sum_{i=1}^{C} (p_i-q_i)^2 + \lambda(\sum_{i=1}^{C} q_i - 1)\] Điều kiện đạt cực tiểu:

\[\Delta_{q_1, q_2, \dots, q_n,\lambda}\mathcal{L}(q_1, q_2, \dots , q_n, \lambda) = 0 \Leftrightarrow \left\{ \begin{matrix} -2~q_i(p_i-q_i) + \lambda &=& 0 ~~, \forall i = 1,\dots, n \\ \sum_{i=1}^{C} q_i &=& 1 \end{matrix} \right.\]

Đặt \(e_i = p_i - q_i, \forall i\). Ta có \(1 = \sum_{i=1}^{C} q_i = \sum_{i=1}^{C} p_i \Rightarrow \sum_{i=1}^C e_i = 0\). Từ phương trình thứ nhất của hệ phương trình suy ra \(q_ie_i = \frac{\lambda}{2}\). Giả sử tồn tại \(e_i \neq 0\), do \(\sum_{i=1}^C e_i = 0\) nên luôn tồn tại một cặp \((e_m, e_n), m \neq n\) sao cho \[\left\{ \begin{matrix} e_m > 0\\ e_n < 0 \end{matrix}\right. \Rightarrow q_me_mq_ne_n = \frac{\lambda^2}{4}< 0\]

Điều này là vô lý cho nên ta suy ra \(e_i = p_i - q_i = 0 \Leftrightarrow p_i = q_i, \forall i\) hay nghiệm của phương trình entrophy và norm bậc 2 là đồng nhất.

Tuy nhiên trong các bài toán xác xuất chúng ta thường sử dụng hàm entrophy hơn vì entrophy có sự chênh lệch giữa giá trị hàm mất mát tại điểm cực tiểu và tại một điểm khác là lớn hơn nên tốc độ hội tụ của gradient descent là nhanh hơn.

Đối với trường hợp có 2 class thì tại data point thứ i xác xuất thực tế để rơi vào class 1 là \(y_{i}\) và xác xuất ước lượng từ hàm sigmoid là \(z_i\). Còn lại sẽ là xác xuất thực tế để rơi vào class 0 là \(1-y_{i}\) và xác xuất ước lượng là \(1-z_i\) Khi đó hàm entropy sẽ có dạng:

\[H(\mathbf{y},\mathbf{z}) = -[\mathbf{y}\log(\mathbf{z}) + \mathbf{(1-y)}\log(\mathbf{(1-z)})] = - \sum_{i=1}^N [y_{i} \log(z_i)+ (1-y_{i}) \log(1-z_i)]\]

Ta nhận thấy đây chính là hàm mất mát \(J(\mathbf{w})\) bài toán logistic regression. Như vậy hàm mất mát logistic regression là trường hợp đặc biệt của entrophy với số class C = 2 hay còn gọi là binary classification.

4 Xây dựng thuật toán trên R

Các hàm số cơ bản:
sigmoid: Hàm sigmoid trả về giá trị ước lượng của xác xuất.
norm2: Tính norm bậc 2 của vector.
logistic_regression: Hồi qui logistic dựa trên gradient descent.

X <- matrix(c(rep(1,length(ExpYears)),ExpYears), ncol = length(ExpYears), byrow = TRUE)
y <- Boss

sigmoid <- function(x){
  1/(1+exp(-x))
}

norm2 <- function(x){
  sum(x^2)
}

logistic_regression <- function(X, y, w_init, eta, tot = 1e-4, max_count = 10000){
  w = list(w_init)
  N = as.numeric(ncol(X))
  d = as.numeric(nrow(X))
  count = 0
  check_after = 20
  while(count <= max_count){
    #mix_data
    data = sample(c(1:N),N)
    #Cập nhật w theo từng điểm dữ liệu
    for (i in data) {
      xi = X[,i] 
      yi = y[i]
      zi = sigmoid(w[[length(w)]] %*% xi %>% as.numeric())
      w_new = w[[length(w)]] + eta*(yi-zi)*xi
      #Kiểm tra điều kiện để break vòng lặp
      if(count %% check_after == 0 & norm2(w_new-w[[length(w)]]) < tot){
        break()
      }
      w[[length(w)+1]] <- w_new
    }
    count = count + 1
  }
  list(w = w,count = count-1)
}

#Khởi tạo vector hệ số w_init
w_init <- runif(2,0,1)
result <- logistic_regression(X, y, w_init, eta = 0.05)
#Số vòng lặp để dừng
result$count
## [1] 10000
w <- result$w
#Hệ số cuối cùng
w[[length(w)]]
## [1] -3.169615  1.036674

Sử dụng hệ số ước lượng cuối cùng của vòng lặp để dự báo xác xuất trở thành xếp tại mỗi data point.

w[[length(w)]] %*% X %>% apply(2,sigmoid)
##  [1] 0.1059361 0.2504401 0.3594092 0.9752154 0.4336617 0.7265289 0.8822339
##  [8] 0.1808001 0.7265289 0.4336617 0.9834894 0.9940818 0.1659504 0.6127164
## [15] 0.8168903 0.6606248 0.9978932 0.1272377 0.9263613 0.9548011

Biểu diễn hàm hàm sigmoid trên đồ thị:

x <- seq(1,9, by = 0.1)
w_0 <- w[[length(w)]][1]
w_1 <- w[[length(w)]][2]
y <- sigmoid(w_0+w_1*x)
df_sigmoid <- data.frame(x,y)
cutpoint <- -w_0/w_1


df %>% ggplot(aes(x = ExpYears, y = Boss)) +
  geom_point(aes(group = Boss, color = Boss, shape = Boss %>% as.character()), size = 2, show.legend = FALSE) +
  geom_line(aes(x = x, y = y), data = df_sigmoid) +
  geom_point(aes(x = cutpoint, y = 0.5), shape = 2, size = 4) +
  geom_vline(aes(xintercept = cutpoint)) +
  scale_x_continuous(breaks = seq(1,10,1))

Ta có thể thấy đồ thị của phương trình sigmoid thu được có dạng một đường cong có xác xuất tăng dần theo năm kinh nghiệm ExpYears. Khi lựa chọn mốc phân loại là điểm có xác xuất là 0.5 thì điểm này rơi vào vị trí giao nhau giữa đường thẳng y = 0.5 và đường sigmoid được đánh dấu tam giác. Đường boundary giữa 2 class chính là 1 đường thẳng và đi qua điểm tam giác nêu trên. Ta cũng dễ dàng chứng minh được boundary khi sử dụng hàm sigmoid là đường thẳng. Giả sử lựa chọn ngưỡng phân chia các class là 0.5. Khi đó class của một điểm phụ thuộc vào giá trị của hàm sigmoid. Nếu điểm rơi vào class 1 thì:

\[\begin{array}&sigmoid(\mathbf{x_i}) &=& \frac{1}{1+e^{-\mathbf{w}^{T}\mathbf{x_i}}} \geq 0.5 \\ \Leftrightarrow e^{-\mathbf{w}^{T}\mathbf{x_i}} &\leq& 1 \\ \Leftrightarrow \mathbf{w}^{T}\mathbf{x_i} &\geq& 0 ~~ (4) \end{array}\]

Phương trình (4) cho thấy boundary để xác định class 1 là một đường thẳng. Tương tự với trường hợp rơi vào class 0.

Mặc dù có boundary là một đường thẳng như phương hàm sigmoid vẫn có lợi thế hơn thuật toán PLA đó là ước lượng được xác xuất xảy ra của một biến cố và sigmoid vẫn hoạt động trong trường hợp 2 class overlap nhau dẫn tới không tồn tại một boundary phân biệt hoàn toàn 2 nhóm trên thực tế. Vì những lý do này mà trong bài toán phân lớp binary class thì logistic regression thường được sử dụng.