機器學習及演算法
林嶔 (Lin, Chin)
Lesson 6 人工智慧基礎3（MCMC與梯度下降法）
第一節：馬可夫鏈蒙特卡羅法(1)
- 在上節課中，我們提到了函數「mle」在求最大概似估計值時並非使用微分法，那他是用什麼方法呢?
– 我們現在已經了解到一件事情，無論函數型態寫成什麼樣子，最大概似估計法重點就是能讓我們定義出一個求解函數，所以我們把問題簡化一點：我們希望有一個方法，能在某個不特定函數中找出該函數的極值。
- 馬可夫鏈蒙特卡羅法(MCMC)的概念是在函數空間內隨機的給予一個起始點，並且命令隨機移動，若新的一點比舊的一點「樣本機率」來的更高，則直接移動，反之則退回原位
– 由於MCMC是完全隨機的移動，在他移動進入極值時可以利用這個隨機移動的特性去找尋該區域的分布特性
- 我們在這裡介紹MCMC中的經典算法：Metropolis–Hastings algorism
第一節：馬可夫鏈蒙特卡羅法(2)
- 在這裡開始想像一個問題，假定有個方程式\(y = x^2\)，我們現在想要找出x等於多少時y有極值存在
original.fun <- function(x) {return(x^2)}
random.walk <- function(x, val = 0.1) {x + runif(1, min = -val, max = val)}
set.seed(0)
start.value = 5
num.iteration = 1000
x = rep(NA, num.iteration)
for (i in 1:100) {
if (i == 1) {
x[i] <- start.value
} else {
old.x <- x[i-1]
new.x <- random.walk(old.x)
if (original.fun(old.x) < original.fun(new.x)) {
x[i] <- old.x
} else {
x[i] <- new.x
}
}
}
x[71:100]
## [1] 3.482776 3.482776 3.450591 3.450591 3.419928 3.386683 3.381953 3.381953
## [9] 3.381953 3.359951 3.359951 3.359951 3.346883 3.346883 3.326882 3.291952
## [17] 3.291952 3.232490 3.232490 3.156829 3.105927 3.034587 2.982513 2.894300
## [25] 2.894300 2.894300 2.894300 2.894300 2.885355 2.867372
- 雖然真的往最小值移動了，但感覺得出來收斂的速度很慢，這是因為方向隨機造成的困擾
for (i in 101:1000) {
old.x = x[i-1]
new.x = random.walk(old.x)
if (original.fun(old.x) < original.fun(new.x)) {
x[i] = old.x
} else {
x[i] = new.x
}
}
x[971:1000]
## [1] -0.0001650651 -0.0001650651 -0.0001650651 -0.0001650651 -0.0001650651
## [6] -0.0001650651 -0.0001650651 -0.0001650651 -0.0001650651 -0.0001650651
## [11] -0.0001650651 -0.0001650651 -0.0001650651 -0.0001650651 -0.0001650651
## [16] -0.0001650651 -0.0001650651 -0.0001650651 -0.0001650651 -0.0001650651
## [21] -0.0001650651 -0.0001650651 -0.0001650651 -0.0001650651 -0.0001650651
## [26] -0.0001650651 -0.0001650651 -0.0001650651 -0.0001650651 -0.0001650651
- 小到一定程度後就幾乎動不了了，該如何搜尋他的機率分布呢？
第一節：馬可夫鏈蒙特卡羅法(3)
- 我們把它應用來解線性迴歸，但現在我們要用最大概似估計法來做最佳化，所以我們求解的目標改變為求樣本機率最大的地方
– 請至這裡下載範例資料
dat <- read.csv("ECG_train.csv", header = TRUE, fileEncoding = 'CP950', stringsAsFactors = FALSE, na.strings = "")
subdat <- dat[!(dat[,"K"] %in% NA) & !(dat[,"PR"] %in% NA),]
x <- subdat[,"PR"]
y <- subdat[,"K"]
prop.fun <- function(b0, b1, x = x, y = y) {
y.hat <- b0 + b1 * x
res <- y.hat - y
sd.res <- sd(res)
log_p <- dnorm(res, mean = 0, sd = sd.res, log = TRUE)
return(sum(log_p))
}
第一節：馬可夫鏈蒙特卡羅法(4)
- 讓我們為他增加一個巧思，在他移動時如果新的機率比舊的機率大，那就照常移動，但若新的機率比舊的機率小，我們仍然給他一個很小的機率移動，但這個移動的機率取決於新機率與舊機率的差異
start.b0 <- 0
start.b1 <- 0
num.iteration <- 10000
b0.seq <- rep(NA, num.iteration)
b1.seq <- rep(NA, num.iteration)
for (i in 1:num.iteration) {
if (i == 1) {
b0.seq[i] <- start.b0
b1.seq[i] <- start.b1
} else {
b0.seq[i] <- random.walk(b0.seq[i-1], val = 1)
b1.seq[i] <- random.walk(b1.seq[i-1], val = 0.1)
old.log_p <- prop.fun(b0 = b0.seq[i-1], b1 = b1.seq[i-1], x = x, y = y)
new.log_p <- prop.fun(b0 = b0.seq[i], b1 = b1.seq[i], x = x, y = y)
diff.p <- exp(new.log_p - old.log_p)
if (diff.p < runif(1, min = 0, max = 1)) {
b0.seq[i] <- b0.seq[i-1]
b1.seq[i] <- b1.seq[i-1]
}
}
}
第一節：馬可夫鏈蒙特卡羅法(5)
par(mfcol = c(1, 2))
plot(1:num.iteration, b0.seq, type = "l")
plot(1:num.iteration, b1.seq, type = "l")
- 看起來是經過一陣子之後就開始收斂了，並且在答案附近上下震動
– 我們設置一個Burn-In time，看看在開始上下震動之後的結果
burn_in <- 5000
par(mfcol = c(2, 2))
hist(b0.seq[(burn_in+1):num.iteration])
abline(v = mean(b0.seq[(burn_in+1):num.iteration]), col = "blue")
abline(v = 5, col = "red")
plot((burn_in+1):num.iteration, b0.seq[(burn_in+1):num.iteration], type = "l")
abline(h = mean(b0.seq[(burn_in+1):num.iteration]), col = "blue")
abline(h = 5, col = "red")
hist(b1.seq[(burn_in+1):num.iteration])
abline(v = mean(b1.seq[(burn_in+1):num.iteration]), col = "blue")
abline(v = 3, col = "red")
plot((burn_in+1):num.iteration, b1.seq[(burn_in+1):num.iteration], type = "l")
abline(h = mean(b1.seq[(burn_in+1):num.iteration]), col = "blue")
abline(h = 3, col = "red")
第一節：馬可夫鏈蒙特卡羅法(6)
- 雖然沒有非常準，但我們已經求出了他的分布了，我們來看看我們的結果
fit <- lm(y~x)
summary(fit)
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2563 -0.8290 -0.6899 0.4701 6.3229
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.5733390 0.1067549 33.472 <2e-16 ***
## x 0.0008796 0.0006399 1.375 0.169
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.264 on 2820 degrees of freedom
## Multiple R-squared: 0.0006697, Adjusted R-squared: 0.0003153
## F-statistic: 1.89 on 1 and 2820 DF, p-value: 0.1693
mean(b0.seq[(burn_in+1):num.iteration])
## [1] 3.549336
mean(b1.seq[(burn_in+1):num.iteration])
## [1] 0.001074259
sd(b0.seq[(burn_in+1):num.iteration])
## [1] 0.09483694
sd(b1.seq[(burn_in+1):num.iteration])
## [1] 0.000531151
- 在預測結果是不是挺相似的，而這個方法最有趣的地方在於他能取得標準誤！
練習1：修正隨機參數
- 透過馬可夫鏈蒙特卡羅法的演示，你應該能總結許多關鍵了。
– 這是我們第一次完全沒有透過數學方式解「函數極值」，而同樣也能獲得「近似解」
– MCMC似乎沒有什麼太大的限制，我們只需要指定一個預測公式以及優化目標，之後就能輕鬆解出答案了
- 剛剛的分析有發現，其中「random.walk()」中的【val】非常關鍵，我們希望最開始稍微大一點，之後數值不太變化時再小一點，你能想到好方法改進嘛？
練習1答案(1)
- 我們可以透過迴圈參數\(i\)來決定隨機大小，隨著時間的進展，我們會希望隨機越來越小：
start.b0 <- 0
start.b1 <- 0
num.iteration <- 10000
b0.seq <- rep(NA, num.iteration)
b1.seq <- rep(NA, num.iteration)
for (i in 1:num.iteration) {
decay_coef <- sqrt(i / 10 + 1)
if (i == 1) {
b0.seq[i] <- start.b0
b1.seq[i] <- start.b1
} else {
b0.seq[i] <- random.walk(b0.seq[i-1], val = 10 / decay_coef)
b1.seq[i] <- random.walk(b1.seq[i-1], val = 1 / decay_coef)
old.log_p <- prop.fun(b0 = b0.seq[i-1], b1 = b1.seq[i-1], x = x, y = y)
new.log_p <- prop.fun(b0 = b0.seq[i], b1 = b1.seq[i], x = x, y = y)
diff.p <- exp(new.log_p - old.log_p)
if (diff.p < runif(1, min = 0, max = 1)) {
b0.seq[i] <- b0.seq[i-1]
b1.seq[i] <- b1.seq[i-1]
}
}
}
練習1答案(2)
burn_in <- 5000
par(mfcol = c(2, 2))
hist(b0.seq[(burn_in+1):num.iteration])
abline(v = mean(b0.seq[(burn_in+1):num.iteration]), col = "blue")
abline(v = 5, col = "red")
plot((burn_in+1):num.iteration, b0.seq[(burn_in+1):num.iteration], type = "l")
abline(h = mean(b0.seq[(burn_in+1):num.iteration]), col = "blue")
abline(h = 5, col = "red")
hist(b1.seq[(burn_in+1):num.iteration])
abline(v = mean(b1.seq[(burn_in+1):num.iteration]), col = "blue")
abline(v = 3, col = "red")
plot((burn_in+1):num.iteration, b1.seq[(burn_in+1):num.iteration], type = "l")
abline(h = mean(b1.seq[(burn_in+1):num.iteration]), col = "blue")
abline(h = 3, col = "red")
mean(b0.seq[(burn_in+1):num.iteration])
## [1] 3.472149
mean(b1.seq[(burn_in+1):num.iteration])
## [1] 0.001550942
sd(b0.seq[(burn_in+1):num.iteration])
## [1] 0.1795203
sd(b1.seq[(burn_in+1):num.iteration])
## [1] 0.001078361
第二節：梯度下降法(1)
- MCMC的確是可以拿來做優化問題的解法，但你似乎看出了問題，那就是速度很慢
– 這個問題是因為隨機移動造成的，如果想要解決這個問題，那我們勢必需要有方法導引他往正確的方向走。
- 在這裡我們隆重介紹『梯度下降法』。首先我們要先定義何謂『梯度』?所謂的『梯度』其實就是『斜率』（注意，這個定義並不精確，但為了省略太多複雜的數學語言，我們暫且使用這個定義）。在這個定義之下，『梯度下降法』意思就是我們在『求解極值』的過程中，隨著『梯度』進行移動，從而找到極值的過程。
– 還是很難理解吧，我們來想想這個函數的求解過程：
\[f(x) = x^{2} + 2x + 1\]
– 我們先把他的導函數寫下來：
\[\frac{\partial}{\partial x} f(x) = 2x + 2\] – 他的意思是說在任何一個點的切線斜率是\(2x + 2\)，而斜率的意思就是說「x每增加一個單位，y所改變的量」
第二節：梯度下降法(2)
- 既然都要求導函數了，為什麼不直接像第五課那樣找導函數等於0的位置硬解呢?
\[\frac{\partial}{\partial x} f(x) = 2x + 2 = 0\]
\[x = -1\]
為什麼我們能夠利用『微分』求函數的極值?這邊大家可能要複習一下基本觀念，對一個『函數』進行『微分』所獲得的『導函數』其實就是該函數的『切線斜率函數』，而『切線斜率函數』等於0的位置就暗示著函數不經過一系列的上升/下降後停止變化，那當然這個位置就是極值所在。
然而，剛剛的求極值過程中有一個非常討厭的部分，那就是要求一個「一元一次方程式」，而當函數複雜一點，我們將要求「N元M次聯立方程式」的答案，那將會讓整個過程異常複雜，所以我們要尋求其他解決方案。
第二節：梯度下降法(3)
\[x_{\left(epoch:0\right)} = 10\]
- 接著，我們定義一下梯度下降法的公式(\(lr\)為學習率，一般我們會給一個很小的值，如下面的範例我們將使用0.05)：
\[x_{\left(epoch:t\right)} = x_{\left(epoch:t - 1\right)} - lr \cdot \frac{\partial}{\partial x}f(x_{\left(epoch:t - 1\right)})\] - 由於剛剛函數的導函數為「\(2x + 2\)」，我們可以將式子帶入運算：
\[
\begin{align}
x_{\left(epoch:1\right)} & = x_{\left(epoch:0\right)} - lr \cdot \frac{\partial}{\partial x}f(x_{\left(epoch:0\right)}) \\
& = 10 - lr \cdot \frac{\partial}{\partial x}f(10) \\
& = 10 - 0.05 \cdot (2\cdot10+2)\\
& = 8.9
\end{align}
\]
第二節：梯度下降法(4)
\[
\begin{align}
x_{\left(epoch:2\right)} & = x_{\left(epoch:1\right)} - lr \cdot \frac{\partial}{\partial x}f(x_{\left(epoch:1\right)}) \\
& = 8.9 - lr \cdot \frac{\partial}{\partial x}f(8.9) \\
& = 8.9 - 0.05 \cdot (2\cdot8.9+2)\\
& = 7.91
\end{align}
\]
\[
\begin{align}
x_{\left(epoch:3\right)} & = 7.91 - 0.891 = 7.019 \\
x_{\left(epoch:4\right)} & = 7.019 - 0.8019 = 6.2171 \\
x_{\left(epoch:5\right)} & = 6.2171 - 0.72171 = 5.49539 \\
& \dots \\
x_{\left(epoch:\infty\right)} & = -1
\end{align}
\]
第二節：梯度下降法(5)
- 讓我們用R語言實現一下梯度下降法，我們稍微簡化一下問題，求這個函數的極值：
\[f(x) = x^{2}\]
original.fun = function(x) {
return(x^2)
}
differential.fun = function(x) {
return(2*x)
}
start.value = 5
learning.rate = 0.1
num.iteration = 1000
result.x = rep(NA, num.iteration)
for (i in 1:num.iteration) {
if (i == 1) {
result.x[1] = start.value
} else {
result.x[i] = result.x[i-1] - learning.rate * differential.fun(result.x[i-1])
}
}
print(tail(result.x, 1))
[1] 7.68895e-97
第二節：梯度下降法(6)
- 我們試著把\(lr\)調高，也就是剛剛程式碼中的learning.rate，看看會發生甚麼事情：
- 看來learning.rate是一個重要的參數，設的太大可能無法收斂，但設的太小幅度又太慢。
– 在使用梯度下降法時，原則上learning.rate不宜設置太大，但可以觀察收斂速度，若收斂速度太慢再適當的調整為佳。
第二節：梯度下降法(7)
- 真實世界的問題往往沒有這麼簡單，現在我們試著來解一個更複雜的極值問題，他一次牽涉到了2個元素：
\[f(a, b) = 4a^{2}- 4a + b^{2} - 2b + 5\]
\[
\begin{align}
f(a, b) & = (4a^{2} - 4a + 1) + (b^{2} - 2b + 1) + 3 \\
& = (2a - 1)^{2} + (b - 1)^{2} + 3
\end{align}
\]
- 好，這下我們知道答案了，當\(a = 0.5\)並且\(b = 1\)時，此函數會有最小值\(3\)。
第二節：梯度下降法(8)
- 現在讓我們試著用梯度下降法求解，這裡我們要使用偏微分求導函數：
\[
\begin{align}
\frac{\partial}{\partial a} f(a, b) & = 8a - 4 \\
\frac{\partial}{\partial b} f(a, b) & = 2b - 2
\end{align}
\]
\[
\begin{align}
a_{\left(epoch:t\right)} & = a_{\left(epoch:t - 1\right)} - lr \cdot \frac{\partial}{\partial a}f(a_{\left(epoch:t - 1\right)}, b_{\left(epoch:t - 1\right)}) \\
b_{\left(epoch:t\right)} & = b_{\left(epoch:t - 1\right)} - lr \cdot \frac{\partial}{\partial b}f(a_{\left(epoch:t - 1\right)}, b_{\left(epoch:t - 1\right)})
\end{align}
\]
- 給個隨機起始值，然後開始…(\(lr\)設為0.1)
\[
\begin{align}
a_{\left(epoch:0\right)} & = 0 \\
b_{\left(epoch:0\right)} & = 0 \\\\
a_{\left(epoch:1\right)} & = a_{\left(epoch:0\right)} - lr \cdot \frac{\partial}{\partial a}f(a_{\left(epoch:0\right)}, b_{\left(epoch:0\right)}) \\
& = 0 - lr \cdot \frac{\partial}{\partial a}f(0, 0) \\
& = 0 - 0.1 \cdot (8\cdot0-4)\\
& = 0.4 \\
b_{\left(epoch:1\right)} & = b_{\left(epoch:0\right)} - lr \cdot \frac{\partial}{\partial b}f(a_{\left(epoch:0\right)}, b_{\left(epoch:0\right)}) \\
& = 0 - lr \cdot \frac{\partial}{\partial b}f(0, 0) \\
& = 0 - 0.1 \cdot (2\cdot0-2)\\
& = 0.2 \\
\end{align}
\]
第二節：梯度下降法(9)
- 讓我們重複這個過程幾次，再度體驗梯度下降法的威力…
\[
\begin{align}
a_{\left(epoch:2\right)} & = 0.4 + 0.08 = 0.48\\
b_{\left(epoch:2\right)} & = 0.2 + 0.16 = 0.36\\\\
a_{\left(epoch:3\right)} & = 0.48 + 0.016 = 0.496\\
b_{\left(epoch:3\right)} & = 0.36 + 0.128 = 0.488\\\\
a_{\left(epoch:4\right)} & = 0.496 + 0.0032 = 0.4992\\
b_{\left(epoch:4\right)} & = 0.488 + 0.1024 = 0.5904\\\\
& \dots \\\\
a_{\left(epoch:\infty\right)} & = 0.5\\
b_{\left(epoch:\infty\right)} & = 1
\end{align}
\]
練習2：利用梯度下降法求解線性迴歸
- 我們試著來用梯度下降法求解線性迴歸吧，數學的部分讓我幫你們解好。首先再次定義其預測函數及損失函數：
– 預測函數
\[\hat{y_{i}} = f(x) = b_{0} + b_{1}x_{i}\]
– 損失函數
\[loss = diff(y, \hat{y}) = \frac{{1}}{2n}\sum \limits_{i=1}^{n} \left(y_{i} - \hat{y_{i}}\right)^{2}\]
– 整合預測函數與損失函數
\[loss = diff(y, f(x)) = \frac{{1}}{2n}\sum \limits_{i=1}^{n} \left(y_{i} - \left(b_{0} + b_{1}x_{i}\right)\right)^{2}\]
- 我們的求解目標為找出一組特定的\(b_0\)以及\(b_1\)使\(loss\)最小化，因此我們將對\(loss\)做偏微分(過程可以在第4課找到)：
– \(b_0\)以及\(b_1\)的偏導函數
\[
\begin{align}
\frac{\partial}{\partial b_0}loss & = \frac{1}{n} \sum \limits_{i=1}^{n}\left( \hat{y_{i}} - y_{i} \right) \\
\frac{\partial}{\partial b_1}loss & = \frac{1}{n} \sum \limits_{i=1}^{n}\left( \hat{y_{i}} - y_{i} \right) \cdot x_{i}
\end{align}
\]
\[
\begin{align}
b_{0\left(epoch:t\right)} & = b_{0\left(epoch:t - 1\right)} - lr \cdot \frac{\partial}{\partial b_{0}}loss \\
b_{1\left(epoch:t\right)} & = b_{1\left(epoch:t - 1\right)} - lr \cdot \frac{\partial}{\partial b_{1}}loss
\end{align}
\]
x <- c(1, 2, 3, 4, 5)
y <- c(6, 7, 9, 8, 10)
original.fun <- function(b0, b1, x = x, y = y) {
y.hat = b0 + b1 * x
return(sum((y.hat - y)^2)/2/length(x))
}
differential.fun.b0 <- function(b0, b1, x = x, y = y) {
y.hat = b0 + b1 * x
return(sum(y.hat - y)/length(x))
}
differential.fun.b1 <- function(b0, b1, x = x, y = y) {
y.hat = b0 + b1 * x
return(sum((y.hat - y)*x)/length(x))
}
model = lm(y~x)
print(model)
##
## Call:
## lm(formula = y ~ x)
##
## Coefficients:
## (Intercept) x
## 5.3 0.9
練習2答案
lr <- 0.1
num.iteration <- 1000
ans_b0 = rep(0, num.iteration)
ans_b1 = rep(0, num.iteration)
for (i in 2:num.iteration) {
ans_b0[i+1] <- ans_b0[i] - lr * differential.fun.b0(b0 = ans_b0[i], b1 = ans_b1[i], x = x, y = y)
ans_b1[i+1] <- ans_b1[i] - lr * differential.fun.b1(b0 = ans_b0[i], b1 = ans_b1[i], x = x, y = y)
}
print(tail(ans_b0, 1))
[1] 5.3
[1] 0.9000001
課程小結
梯度下降法是目前最有效率的求解方式，他的應用也最為廣泛，限制只有一個，那就是要保證每個參數都是有導函數的！
在人工智慧課程開始至今，我們已經不知不覺熟悉的一個人工智慧模型的訓練過程，我們再重新整理一下。
所謂的人工智慧，說白了就是個數學方程式，在數學上我們可以假定Input是\(x\)物件，而Output是\(y\)物件，那我們可以建構一個『預測函數』\(f()\)來對其進行預測，而這個預測函數是可以很複雜的。
\[\hat{y} = f(x)\] - 但凡預測就會存在誤差，我們可以再定義一個『損失函數』\(diff()\)，而這個損失函數可以由最大概似估計法決定：
\[loss = diff(y, \hat{y})\] - 接著就是使用梯度下降法求解，最終我們就能求得一組參數，讓\(loss\)最小化，從而得到精準的預測。
- 現在的深度神經網路完全就是基於這項技術而生的，而這整個流程我們還會用來解非常多方程式，大家就藉著這整套邏輯開始慢慢探索資料科學的世界吧！