2019/9/15

以 python 實作迴歸與分類程式

迴歸


學習資料


取得學習資料文字檔 click.csv


x,y
235,591
216,539
148,413
35,310
85,308
...

先利用 matplotlib 繪製到圖表上


import numpy as np
import matplotlib.pyplot as plt

train = np.loadtxt('click.csv', delimiter=',', skiprows=1)
train_x = train[:,0]
train_y = train[:,1]

plt.plot(train_x, train_y, 'o')
plt.show()


另外針對原始的學習資料,進行標準化(z-score正規化),也就是將資料平均轉換為 0,分散轉換為1。其中 𝜇 是所有資料的平均,𝜎 是所有資料的標準差。這樣處理後,會讓參數收斂更快。


\(z^{(i)} = \frac{x^{(i)} - 𝜇}{𝜎}\)


import numpy as np
import matplotlib.pyplot as plt

train = np.loadtxt('click.csv', delimiter=',', skiprows=1)
train_x = train[:,0]
train_y = train[:,1]

# 標準化
mu = train_x.mean()
sigma = train_x.std()
def standardize(x):
    return (x - mu) / sigma

train_z = standardize(train_x)

plt.plot(train_z, train_y, 'o')
plt.show()


一次函數


先使用一次目標函數 \(f_𝜃(x)\)


\({f_𝜃(x)=𝜃_0+𝜃_1x}​\)


\({E(𝜃)= \frac{1}{2} \sum_{i=1}^{n}( y^{(i)} - f_𝜃(x^{(i)})^2 }​\)


\(𝜃_0, 𝜃_1​\) 可任意選擇初始值


\(𝜃_0, 𝜃_1\) 的參數更新式為


\(𝜃_0 := 𝜃_0 - 𝜂 \sum_{i=1}^{n}( f_𝜃(x^{(i)} )-y^{(i)} )\)


\(𝜃_1 := 𝜃_1 - 𝜂 \sum_{i=1}^{n}( f_𝜃(x^{(i)} )-y^{(i)} )x^{(i)}\)


用這個方法,就可以找出正確的 \(𝜃_0, 𝜃_1\)


其中 𝜂 是任意數值,先設定為 \(10^{-3}\) 試試看。一般來說,會指定要處理的次數,有時會比較參數更新前後,目標函數的值,如果差異不大,就直接結束。另外 \(𝜃_0, 𝜃_1\) 必須同時一起更新。


import numpy as np
import matplotlib.pyplot as plt

# 載入學習資料
train = np.loadtxt('click.csv', delimiter=',', dtype='int', skiprows=1)
train_x = train[:,0]
train_y = train[:,1]

# 標準化
mu = train_x.mean()
sigma = train_x.std()
def standardize(x):
    return (x - mu) / sigma

train_z = standardize(train_x)

# 任意選擇初始值
theta0 = np.random.rand()
theta1 = np.random.rand()

# 預測函數
def f(x):
    return theta0 + theta1 * x

# 目標函數 E(𝜃)
def E(x, y):
    return 0.5 * np.sum((y - f(x)) ** 2)

# 學習率
ETA = 1e-3

# 誤差
diff = 1

# 更新次數
count = 0

# 重複學習
error = E(train_z, train_y)
while diff > 1e-2:
    # 暫存更新結果
    tmp_theta0 = theta0 - ETA * np.sum((f(train_z) - train_y))
    tmp_theta1 = theta1 - ETA * np.sum((f(train_z) - train_y) * train_z)

    # 更新參數
    theta0 = tmp_theta0
    theta1 = tmp_theta1

    # 計算誤差
    current_error = E(train_z, train_y)
    diff = error - current_error
    error = current_error

    # log
    count += 1
    log = '{}次數: theta0 = {:.3f}, theta1 = {:.3f}, 誤差 = {:.4f}'
    print(log.format(count, theta0, theta1, diff))

# 繪製學習資料與預測函數的直線
x = np.linspace(-3, 3, 100)
plt.plot(train_z, train_y, 'o')
plt.plot(x, f(x))
plt.show()

測試結果


391次數: theta0 = 428.991, theta1 = 93.444, 誤差 = 0.0109
392次數: theta0 = 428.994, theta1 = 93.445, 誤差 = 0.0105
393次數: theta0 = 428.997, theta1 = 93.446, 誤差 = 0.0101
394次數: theta0 = 429.000, theta1 = 93.446, 誤差 = 0.0097


驗證


可輸入 x 預測點擊數,但因為剛剛有將學習資料正規化,預測資料也必須正規化


>>> f(standardize(100))
370.96741051658194
>>> f(standardize(500))
928.9775823086377

二次多項式迴歸


\(f_𝜃(x) = 𝜃_0 + 𝜃_1x + 𝜃_2x^2\) 要增加 \( 𝜃_2\) 這個參數


目標的誤差函數 \({E(𝜃)= \frac{1}{2} \sum_{i=1}^{n}( y^{(i)} - f_𝜃(x^{(i)})^2 }​\)


因為有多筆學習資料,可將資料以矩陣方式處理


\( X = \begin{bmatrix}
(x^{(1)})^T\\
(x^{(2)})^T\\
\cdot \\
\cdot \\
(x^{(n)})^T \\
\end{bmatrix}
= \begin{bmatrix}
1 & x^{(1)} & (x^{(1)})^2 \\
1 & x^{(2)} & (x^{(2)})^2 \\
\cdot \\
\cdot \\
1 & x^{(n)} & (x^{(n)})^2 \\
\end{bmatrix} ​\)


\(f_𝜃(x) = \begin{bmatrix}
1 & x^{(1)} & (x^{(1)})^2 \\
1 & x^{(2)} & (x^{(2)})^2 \\
\cdot \\
\cdot \\
1 & x^{(n)} & (x^{(n)})^2 \\
\end{bmatrix} \begin{bmatrix}
𝜃_0 \\
𝜃_1 \\
𝜃_2 \\
\end{bmatrix}
= \begin{bmatrix}
𝜃_0 + 𝜃_1 x^{(1)} + 𝜃_2 (x^{(1)})^2\\
𝜃_0 + 𝜃_1 x^{(2)} + 𝜃_2 (x^{(2)})^2\\
\cdot \\
\cdot \\
𝜃_0 + 𝜃_1 x^{(n)} + 𝜃_2 (x^{(n)})^2\\
\end{bmatrix}\)


第j 項參數的更新式定義為


\(𝜃_j := 𝜃_j - 𝜂 \sum_{i=1}^{n}( f_𝜃(x^{(i)} )-y^{(i)} )x_j^{(i)}​\)


可將 \( ( f_𝜃(x^{(i)} )-y^{(i)} ) ​\) 以及 \(x_j^{(i)}​\) 這兩部分各自以矩陣方式處理


\( f= \begin{bmatrix}
( f_𝜃(x^{(1)} )-y^{(1)} )\\
( f_𝜃(x^{(2)} )-y^{(2)} )\\
\cdot \\
\cdot \\
( f_𝜃(x^{(n)} )-y^{(n)} ) \\
\end{bmatrix} \)


\( x_0 = \begin{bmatrix}
x_0^{(1)} \\
x_0^{(2)}\\
\cdot \\
\cdot \\
x_0^{(n)} \\
\end{bmatrix} \)


\( \sum_{i=1}^{n}( f_𝜃(x^{(i)} )-y^{(i)} )x_0^{(i)} = f^Tx_0 \)


分別考慮三個參數


\( x_0 = \begin{bmatrix}
x_0^{(1)} \\
x_0^{(2)}\\
\cdot \\
\cdot \\
x_0^{(n)} \\
\end{bmatrix} ,
x_1 = \begin{bmatrix}
x^{(1)} \\
x^{(2)}\\
\cdot \\
\cdot \\
x^{(n)} \\
\end{bmatrix} ,
x_2 = \begin{bmatrix}
(x^{(1)})^2 \\
(x^{(2)})^2\\
\cdot \\
\cdot \\
(x^{(n)})^2 \\
\end{bmatrix}\)


\( X = \begin{bmatrix}
x_0 & x_1 & x_2
\end{bmatrix}
= \begin{bmatrix}
1 & x^{(1)} & (x^{(1)})^2 \\
1 & x^{(2)} & (x^{(2)})^2\\
\cdot \\
\cdot \\
1 & x^{(n)} & (x^{(n)})^2 \\
\end{bmatrix} \)


使用 \( f^TX\) 就可以一次更新三個參數


import numpy as np
import matplotlib.pyplot as plt

# 讀取學習資料
train = np.loadtxt('click.csv', delimiter=',', dtype='int', skiprows=1)
train_x = train[:,0]
train_y = train[:,1]

# 標準化
mu = train_x.mean()
sigma = train_x.std()
def standardize(x):
    return (x - mu) / sigma

train_z = standardize(train_x)

# 任意初始值
theta = np.random.rand(3)

# 學習資料轉換為矩陣
def to_matrix(x):
    return np.vstack([np.ones(x.size), x, x ** 2]).T

X = to_matrix(train_z)

# 預測函數
def f(x):
    return np.dot(x, theta)

# 目標函數
def E(x, y):
    return 0.5 * np.sum((y - f(x)) ** 2)

# 學習率
ETA = 1e-3

# 誤差
diff = 1

# 更新次數
count = 0

# 重複學習
error = E(X, train_y)
while diff > 1e-2:
    # 更新參數
    theta = theta - ETA * np.dot(f(X) - train_y, X)

    # 計算誤差
    current_error = E(X, train_y)
    diff = error - current_error
    error = current_error

    # log
    count += 1
    log = '{}次: theta = {}, 誤差 = {:.4f}'
    print(log.format(count, theta, diff))

# 繪製學習資料與預測函數
x = np.linspace(-3, 3, 100)
plt.plot(train_z, train_y, 'o')
plt.plot(x, f(to_matrix(x)))
plt.show()




也可以將重複停止的條件,改為均方誤差


目標的誤差函數 \({E(𝜃)= \frac{1}{n} \sum_{i=1}^{n}( y^{(i)} - f_𝜃(x^{(i)})^2 }\)


import numpy as np
import matplotlib.pyplot as plt

# 讀取學習資料
train = np.loadtxt('click.csv', delimiter=',', dtype='int', skiprows=1)
train_x = train[:,0]
train_y = train[:,1]

# 標準化
mu = train_x.mean()
sigma = train_x.std()
def standardize(x):
    return (x - mu) / sigma

train_z = standardize(train_x)

# 任意初始值
theta = np.random.rand(3)

# 學習資料轉換為矩陣
def to_matrix(x):
    return np.vstack([np.ones(x.size), x, x ** 2]).T

X = to_matrix(train_z)

# 預測函數
def f(x):
    return np.dot(x, theta)

# 目標函數
def MSE(x, y):
    return ( 1 / x.shape[0] * np.sum( (y-f(x)))**2 )

# 學習率
ETA = 1e-3

# 誤差
diff = 1

# 更新次數
count = 0

# 均方誤差的歷史資料
errors = []

# 重複學習
errors.append( MSE(X, train_y) )
while diff > 1e-2:
    # 更新參數
    theta = theta - ETA * np.dot(f(X) - train_y, X)

    # 計算誤差
    errors.append( MSE(X, train_y) )
    diff = errors[-2] - errors[-1]

    # log
    count += 1
    log = '{}次: theta = {}, 誤差 = {:.4f}'
    print(log.format(count, theta, diff))

# 繪製重複次數 與誤差的關係
x = np.arange(len(errors))
plt.plot(x, errors)
plt.show()


隨機梯度下降法


隨機選擇一項學習資料,套用在參數的更新上,例如選擇第 k 項。


\(𝜃_j := 𝜃_j - 𝜂 ( f_𝜃(x^{(k)} )-y^{(k)} )x_j^{(k)}\)


import numpy as np
import matplotlib.pyplot as plt

# 載入學習資料
train = np.loadtxt('click.csv', delimiter=',', dtype='int', skiprows=1)
train_x = train[:,0]
train_y = train[:,1]

# 標準化
mu = train_x.mean()
sigma = train_x.std()
def standardize(x):
    return (x - mu) / sigma

train_z = standardize(train_x)

# 任意選擇初始值
theta = np.random.rand(3)

# 學習資料轉換為矩陣
def to_matrix(x):
    return np.vstack([np.ones(x.size), x, x ** 2]).T

X = to_matrix(train_z)

# 預測函數
def f(x):
    return np.dot(x, theta)

# 均方差
def MSE(x, y):
    return (1 / x.shape[0]) * np.sum((y - f(x)) ** 2)

# 學習率
ETA = 1e-3

# 誤差
diff = 1

# 更新次數
count = 0

# 重複學習
error = MSE(X, train_y)
while diff > 1e-2:
    # 排列學習資料所需的隨機排列
    p = np.random.permutation(X.shape[0])
    # 將學習資料以隨機方式取出,並用隨機梯度下降法 更新參數
    for x, y in zip(X[p,:], train_y[p]):
        theta = theta - ETA * (f(x) - y) * x

    # 計算跟前一個誤差的差距
    current_error = MSE(X, train_y)
    diff = error - current_error
    error = current_error

    # log
    count += 1
    log = '{}回目: theta = {}, 差分 = {:.4f}'
    print(log.format(count, theta, diff))

# 列印結果
x = np.linspace(-3, 3, 100)
plt.plot(train_z, train_y, 'o')
plt.plot(x, f(to_matrix(x)))
plt.show()

多元迴歸


如果要處理多元迴歸,就跟多項式迴歸一樣改用矩陣,但在多元迴歸中要注意,要對所有變數 \(x_1, x_2, x_3\)都進行標準化。


\(z_1^{(i)} = \frac{x_1^{(i)} - 𝜇_1}{𝜎_1} \)


\(z_2^{(i)} = \frac{x_2^{(i)} - 𝜇_2}{𝜎_2} \)


\(z_3^{(i)} = \frac{x_3^{(i)} - 𝜇_3}{𝜎_3} \)


分類(感知器)


使用 images1.csv 資料


x1,x2,y
153,432,-1
220,262,-1
118,214,-1
474,384,1
485,411,1
233,430,-1
...

先將原始資料標記在圖表上,y=1 用圓圈,y=-1 用


import numpy as np
import matplotlib.pyplot as plt

# 載入學習資料
train = np.loadtxt('images1.csv', delimiter=',', skiprows=1)
train_x = train[:,0:2]
train_y = train[:,2]

# 繪圖
x1 = np.arange(0, 500)
plt.plot(train_x[train_y ==  1, 0], train_x[train_y ==  1, 1], 'o')
plt.plot(train_x[train_y == -1, 0], train_x[train_y == -1, 1], 'x')
plt.savefig('1.png')


  • 識別函數 \(f_w(x)\) 就是給定向量 \(x\) 後,回傳 1 或 -1 的函數,用來判斷橫向或縱向。

\(f_w(x) = \left\{\begin{matrix} 1 \quad (w \cdot x \geq 0) \\ -1 \quad (w \cdot x < 0) \end{matrix}\right.\)


  • 權重更新式

\(w := \left\{\begin{matrix} w + y^{(i)}x^{(i)} \quad (f_w(x) \neq y^{(i)}) \\ w \quad \quad \quad \quad (f_w(x) = y^{(i)}) \end{matrix}\right.\)


感知器使用精度作為停止的標準比較好,但目前先直接設定訓練次數


最後繪製以權重向量為法線的直線方程式


\(w \cdot x = w_1x_1 + w_2x_2 = 0​\)


\(x_2 = - \frac{w_1}{w2} x_1​\)


import numpy as np
import matplotlib.pyplot as plt

# 載入學習資料
train = np.loadtxt('images1.csv', delimiter=',', skiprows=1)
train_x = train[:,0:2]
train_y = train[:,2]

# 任意初始值
w = np.random.rand(2)

# 識別函數,判斷矩形是橫向或縱向
def f(x):
    if np.dot(w, x) >= 0:
        return 1
    else:
        return -1

# 重複次數
epoch = 10

# 更新次數
count = 0

# 學習權重
for _ in range(epoch):
    for x, y in zip(train_x, train_y):
        if f(x) != y:
            w = w + y * x

            # log
            count += 1
            print('{}次數: w = {}'.format(count, w))

# 繪圖
x1 = np.arange(0, 500)
plt.plot(train_x[train_y ==  1, 0], train_x[train_y ==  1, 1], 'o')
plt.plot(train_x[train_y == -1, 0], train_x[train_y == -1, 1], 'x')
plt.plot(x1, -w[0] / w[1] * x1, linestyle='dashed')
plt.savefig("1.png")


驗證


python -i classification1_perceptron.py
>>> f([200,100])
1
>>> f([100,200])
-1

分類(邏輯迴歸)


邏輯迴歸要先修改學習資料,橫向為 1 ,縱向為 0


x1,x2,y
153,432,0
220,262,0
118,214,0
474,384,1
485,411,1
...

預測函數就是 S 函數


\(f_𝜃(x) = \frac{1}{1 + exp(-𝜃^Tx)}\)


參數更新式為


\(𝜃_j := 𝜃_j - 𝜂 \sum_{i=1}^{n}( f_𝜃(x^{(i)}) - y^{(i)} )x_j^{(i)}\)


可用矩陣處理,轉換時要加上 \(x_0\),且設定為 1,如果當 \(f_𝜃(x) \geq 0.5\),也就是 \(𝜃^T x >0​\) ,就判定為橫向。


將 \(Q^Tx = 0 \) 整理後,就可得到一條直線


\(Q^Tx = 𝜃_0x_0 + 𝜃_1x_1 + 𝜃_2x_2 = 𝜃_0 +𝜃_1x_1+𝜃_2x_2 =0\)


\(x_2 = - \frac{𝜃_0 + 𝜃_1x_2}{𝜃_2}​\)


import numpy as np
import matplotlib.pyplot as plt

# 載入學習資料
train = np.loadtxt('images2.csv', delimiter=',', skiprows=1)
train_x = train[:,0:2]
train_y = train[:,2]

# 任意初始值
theta = np.random.rand(3)

# 以平均及標準差進行標準化
mu = train_x.mean(axis=0)
sigma = train_x.std(axis=0)
def standardize(x):
    return (x - mu) / sigma

train_z = standardize(train_x)

# 轉換為矩陣,加上 x0
def to_matrix(x):
    x0 = np.ones([x.shape[0], 1])
    return np.hstack([x0, x])

X = to_matrix(train_z)

# 預測函數 S函數
def f(x):
    return 1 / (1 + np.exp(-np.dot(x, theta)))

# 識別函數
def classify(x):
    return (f(x) >= 0.5).astype(np.int)

# 學習率
ETA = 1e-3

# 重複次數
epoch = 5000

# 更新次數
count = 0

# 重複學習
for _ in range(epoch):
    theta = theta - ETA * np.dot(f(X) - train_y, X)

    # log
    count += 1
    print('{}次數: theta = {}'.format(count, theta))

# 繪製圖形
x0 = np.linspace(-2, 2, 100)
plt.plot(train_z[train_y == 1, 0], train_z[train_y == 1, 1], 'o')
plt.plot(train_z[train_y == 0, 0], train_z[train_y == 0, 1], 'x')
plt.plot(x0, -(theta[0] + theta[1] * x0) / theta[2], linestyle='dashed')
# plt.show()
plt.savefig("機器學習4_coding_.png")


驗證


這樣的意思是 200x100 的矩形有 91.6% 的機率會是橫向


>>> f(to_matrix(standardize([[200,100], [100,200]])))
array([0.91604483, 0.03009514])

可再轉化為 1 與 0


>>> classify(to_matrix(standardize([[200,100], [100,200]])))
array([1, 0])

線性不可分離的分類


學習資料為 data3.csv


x1,x2,y
0.54508775,2.34541183,0
0.32769134,13.43066561,0
4.42748117,14.74150395,0
2.98189041,-1.81818172,1
4.02286274,8.90695686,1
2.26722613,-6.61287392,1
-2.66447221,5.05453871,1
-1.03482441,-1.95643469,1
4.06331548,1.70892541,1
2.89053966,6.07174283,0
2.26929206,10.59789814,0
4.68096051,13.01153161,1
1.27884366,-9.83826738,1
-0.1485496,12.99605136,0
-0.65113893,10.59417745,0
3.69145079,3.25209182,1
-0.63429623,11.6135625,0
0.17589959,5.84139826,0
0.98204409,-9.41271559,1
-0.11094911,6.27900499,0

先將學習資料繪製到圖表上看起來無法用一條直線來分類,增加 \(x_1^2​\) 進行分類



參數變成四個,將 \(Q^Tx = 0 ​\) 整理後,就可得到一條曲線


\(Q^Tx = 𝜃_0x_0 + 𝜃_1x_1 + 𝜃_2x_2 +𝜃_3x_1^2 = 𝜃_0 +𝜃_1x_1+𝜃_2x_2 +𝜃_3x_1^2 =0​\)


\(x_2 = - \frac{𝜃_0 + 𝜃_1x_2 +𝜃_3x_1^2}{𝜃_2}\)


import numpy as np
import matplotlib.pyplot as plt

# 載入學習資料
train = np.loadtxt('data3.csv', delimiter=',', skiprows=1)
train_x = train[:,0:2]
train_y = train[:,2]

# 任意初始值
theta = np.random.rand(4)

# 標準化
mu = train_x.mean(axis=0)
sigma = train_x.std(axis=0)
def standardize(x):
    return (x - mu) / sigma

train_z = standardize(train_x)

# 轉換為矩陣,加上 x0, x3
def to_matrix(x):
    x0 = np.ones([x.shape[0], 1])
    x3 = x[:,0,np.newaxis] ** 2
    return np.hstack([x0, x, x3])

X = to_matrix(train_z)

# 預測函數 S函數
def f(x):
    return 1 / (1 + np.exp(-np.dot(x, theta)))

# 識別函數
def classify(x):
    return (f(x) >= 0.5).astype(np.int)

# 學習率
ETA = 1e-3

# 重複次數
epoch = 5000

# 更新次數
count = 0

# 重複學習
for _ in range(epoch):
    theta = theta - ETA * np.dot(f(X) - train_y, X)

    # log
    count += 1
    print('{}次數: theta = {}'.format(count, theta))

# 繪製圖形
x1 = np.linspace(-2, 2, 100)
x2 = -(theta[0] + theta[1] * x1 + theta[3] * x1 ** 2) / theta[2]
plt.plot(train_z[train_y == 1, 0], train_z[train_y == 1, 1], 'o')
plt.plot(train_z[train_y == 0, 0], train_z[train_y == 0, 1], 'x')
plt.plot(x1, x2, linestyle='dashed')
# plt.show()
plt.savefig("機器學習4_coding_.png")


分類的精度,就是在全部的資料中,能夠被正確分類的 TP與 TN 佔的比例,可表示為


\( Accuracy = \frac{TP + TN}{TP+FP+FN+TN} ​\)


# 精度
accuracies = []

# 重複學習
for _ in range(epoch):
    theta = theta - ETA * np.dot(f(X) - train_y, X)

    # 計算精度
    result = classify(X) == train_y
    accuracy = len(result[result ==True]) / len(result)
    accuracies.append(accuracy)

# 繪製圖形
x = np.arange(len(accuracies))
plt.plot(x, accuracies)
plt.savefig("機器學習4_coding_.png")

計算精度,繪製圖表



隨機梯度下降法


import numpy as np
import matplotlib.pyplot as plt

# 載入學習資料
train = np.loadtxt('data3.csv', delimiter=',', skiprows=1)
train_x = train[:,0:2]
train_y = train[:,2]

# 任意初始值
theta = np.random.rand(4)

# 標準化
mu = train_x.mean(axis=0)
sigma = train_x.std(axis=0)
def standardize(x):
    return (x - mu) / sigma

train_z = standardize(train_x)

# 轉換為矩陣,加上 x0, x3
def to_matrix(x):
    x0 = np.ones([x.shape[0], 1])
    x3 = x[:,0,np.newaxis] ** 2
    return np.hstack([x0, x, x3])

X = to_matrix(train_z)

# 預測函數 S函數
def f(x):
    return 1 / (1 + np.exp(-np.dot(x, theta)))

# 識別函數
def classify(x):
    return (f(x) >= 0.5).astype(np.int)

# 學習率
ETA = 1e-3

# 重複次數
epoch = 5000

# 更新次數
count = 0

# 重複學習
for _ in range(epoch):
    # 以隨機梯度下降法更新參數
    p = np.random.permutation(X.shape[0])
    for x, y in zip(X[p,:], train_y[p]):
        theta = theta - ETA * (f(x) - y) * x

    # log
    count += 1
    print('{}次數: theta = {}'.format(count, theta))

# 繪製圖形
x1 = np.linspace(-2, 2, 100)
x2 = -(theta[0] + theta[1] * x1 + theta[3] * x1 ** 2) / theta[2]
plt.plot(train_z[train_y == 1, 0], train_z[train_y == 1, 1], 'o')
plt.plot(train_z[train_y == 0, 0], train_z[train_y == 0, 1], 'x')
plt.plot(x1, x2, linestyle='dashed')
# plt.show()
plt.savefig("機器學習4_coding_.png")


正規化


首先考慮這樣的函數


\(g(x) = 0.1(x^3 + x^2 + x)\)


產生一些雜訊的學習資料,並繪製圖表



import numpy as np
import matplotlib.pyplot as plt

# 原始真正的函數
def g(x):
    return 0.1 * (x ** 3 + x ** 2 + x)

# 適當地利用原本的函數,加上一些雜訊,產生學習資料
train_x = np.linspace(-2, 2, 8)
train_y = g(train_x) + np.random.randn(train_x.size) * 0.05


plt.clf()
x=np.linspace(-2, 2, 100)
plt.plot(train_x, train_y, 'o')
plt.plot(x, g(x), linestyle='dashed')
plt.ylim(-1,2)
plt.savefig("機器學習4_coding_1.png")


# 標準化
mu = train_x.mean()
sigma = train_x.std()
def standardize(x):
    return (x - mu) / sigma

train_z = standardize(train_x)

# 產生學習資料的矩陣 (10次多項式)
def to_matrix(x):
    return np.vstack([
        np.ones(x.size),
        x,
        x ** 2,
        x ** 3,
        x ** 4,
        x ** 5,
        x ** 6,
        x ** 7,
        x ** 8,
        x ** 9,
        x ** 10
    ]).T

X = to_matrix(train_z)

# 參數使用任意初始值
theta = np.random.randn(X.shape[1])

# 預測函數
def f(x):
    return np.dot(x, theta)

# 目標函數
def E(x, y):
    return 0.5 * np.sum((y - f(x)) ** 2)

# 正規化常數
LAMBDA = 0.5

# 學習率
ETA = 1e-4

# 誤差
diff = 1

# 重複學習
error = E(X, train_y)
while diff > 1e-6:
    theta = theta - ETA * (np.dot(f(X) - train_y, X))

    current_error = E(X, train_y)
    diff = error - current_error
    error = current_error

theta1 = theta

# 加上正規化項
theta = np.random.randn(X.shape[1])
diff = 1
error = E(X, train_y)
while diff > 1e-6:
    # 正規化項,因為偏差項不適用於正規化,所以為 0,當 j>0,正規化項為 𝜆 * 𝜃
    reg_term = LAMBDA * np.hstack([0, theta[1:]])
    # 適用於正規化項,更新參數
    theta = theta - ETA * (np.dot(f(X) - train_y, X) + reg_term)

    current_error = E(X, train_y)
    diff = error - current_error
    error = current_error

theta2 = theta

# 繪製圖表
plt.clf()
plt.plot(train_z, train_y, 'o')
z = standardize(np.linspace(-2, 2, 100))
theta = theta1 # 無正規化的結果,虛線
plt.plot(z, f(to_matrix(z)), linestyle='dashed')
theta = theta2 # 有正規化的結果,實線
plt.plot(z, f(to_matrix(z)))
# plt.show()
plt.savefig("機器學習4_coding_2.png")


References


練好機器學習的基本功 範例下載

2019/9/8

機器學習_評估


確認模型的正確性,針對建立的模型,以評估的方法進行機器學習。


迴歸與分類都是定義預測時所需要的函數 \(f_𝜃(x)\),藉由學習資料來找到函數中的 𝜃。方法是將目標函數進行微分,求得參數更新式。但實際上需要的,是預測函數得到的預測值,例如花多少廣告費可得多少點擊率。


我們需要量測函數 \(f_𝜃(x)\) 的正確性(精確度),但像多元迴歸,無法用圖形表示,就需要將機器學習的模型的精度,以定量的方式表示,然後表現其精確度,這就是模型評估。因為參數是透過學習資料修正而來的,對於原本的學習資料來說,參數是正確的,但對於新的資料就不一定了。


交叉驗證 Cross Validation


將學習資料區分為學習以及測試使用,用測試用的資料評估模型,一般來說,學習用的資料會比較多。


在回歸問題中,函數 \(f_𝜃(x)\) 是透過習資料修正而來的,對於原本的學習資料來說,參數是正確的,但對於測試部分的資料就不一定正確了。


假設測試資料有 n 筆,將測試用的資料,透過模型求得結果,再跟原本的實際值比較得到誤差。以下為均方差 MSE (Mean Square Error)


\( \frac{1}{n} \sum_{i=1}^{n} ( y^{(i)} - f_𝜃(x^{(i)} ) )^2 ​\)


當均方差越小,表示模型的精度很高。


分類問題的驗證


因迴歸問題是連續值,可用誤差進行驗證,分類問題是邏輯迴歸,回到矩形是橫向或縱向的問題上,會有四種狀況。


分類結果 原本是橫向 原本是縱向
橫向 正確 錯誤
縱向 錯誤 正確

可將二元分類轉換為這樣的表格


分類結果 + -
+ Positive True Positive (TP) False Positive (FP)
- Negative False Negative (FN) True Negative (TN)

分類的精度,就是在全部的資料中,能夠被正確分類的 TP與 TN 佔的比例,可表示為


\( Accuracy = \frac{TP + TN}{TP+FP+FN+TN} ​\)


ex: 100 筆測試資料,有 80 筆正確


\( Accuracy = \frac{80}{100} = 0.8\)


精確率與回現率


有時候只用精確度評估分類結果會遇到問題。例如當原始資料有大量資料 95 筆為 Negative,只有一點點資料 5 筆為 Positive,如果將全部測試資料都分類為 False 的模型,Accuracy 為 0.95,精確度很高但實際上這是一個錯誤的預測模型。


因此要導入其他評估的指標。


  • 精確率 Precision: 分類為 Positive 的資料中,實際為 Positive 的資料數的比例。值越高,代表分類錯誤的越少。


    \( Precision = \frac{TP}{TP+FP}\)

  • 回現率 Recall: Positive 資料中,實際上分類為 Positive 的資料數。值越高,代表沒有被遺漏,且被正確分類的比例。


    \( Recall = \frac{TP}{TP+FN} \)


通常精確率與回現率,只要有ㄧ個是高的,另一個就會變低。


舉例來說,


資料 個數
Positive 5
Negative 95
評估結果
True Positive 1
False Positive 2
False Negative 4
True Negative 93
精確度 Accuracy 94%
精確率 Precision \(\frac{1}{1+2} = 0.333\)
回現率 Recall \(\frac{1}{1+4} = 0.2\)

F 值 (Fmeasure)


通常精確率與回現率,只要有ㄧ個是高的,另一個就會變低。但直接將兩個平均,也不是好的指標


模型 精確率 回現率 平均
A 0.6 0.39 0.495
B 0.02 1.0 0.51

B 模型是將全部的資料都分類為 Positive,但因為 Negative 也分類為 Positive,所以精確率很低,實際上,B 不是一個好的模型。


Fmeasure 定義如下,只要 Precision 或 Recall 其中一項變低,就會影響到 Fmeasure


\( Fmeasure = \frac{2}{ \frac{1}{Precision} + \frac{1}{Recall} } = \frac{2 \cdot Precision \cdot Recall}{Precision + Recall}\)


模型 精確率 回現率 平均 Fmeasure
A 0.6 0.39 0.495 0.472
B 0.02 1.0 0.51 0.039

F值有時被稱為 F1 值




F值可再加上權重


\(WeightedFmeasure = \frac{ (1+𝛽)^2 \cdot Precision \cdot Recall }{ 𝛽^2 \cdot Precision + Recall }\)


將權重設定為 1 就是原本的 F 值,也就是 F1值




剛剛都是以 TP 為主,考慮精確率與回現率。如果以 TN 為主


\( Precision = \frac{TN}{TN+FN} \)


\( Recall = \frac{TN}{TN+FP} ​\)


當測試資料 Positive 的部分較少,就用 Positive 的 Precision, Recall 來評估。




交叉驗證中,以 K 等分交叉驗證最常見


  • 將學習資料分為 K 筆
  • K-1 筆作為學習用的資料,1 筆作為測試資料
  • 將學習用資料與測試資料,一邊交換,一邊驗證,重複 K 次交叉驗證
  • 最後計算 K 筆精度平均值,視為最終的精度

例如 4 等分



正規化


過適 Overfitting


只跟學習資料吻合的狀態就是 overfitting。如果迴歸中 \(f_𝜃(x)\) 的次方數過度增加,就會造成 overfitting。分類有一樣的問題。


為了避免 overfitting,有以下對應方式


  • 增加學習資料的數量
  • 將模型簡化為較簡單的形式
  • 正規化

正規化


在迴歸分析中的誤差函數為


\({E(𝜃)= \frac{1}{2} \sum_{i=1}^{n}( y^{(i)} - f_𝜃(x^{(i)})^2 }​\)


對該目標函數,再增加正規化的項目 \(R(𝜃) = \frac{𝜆}{2} \sum_{j=1}^{m} 𝜃_j^2​\)


\({E(𝜃)= \frac{1}{2} \sum_{i=1}^{n}( y^{(i)} - f_𝜃(x^{(i)})^2 } + R(𝜃)​\)


對新的目標函數進行最小化,就是正規化


m 是參數的個數,通常對於 \(𝜃_0\) 來說,無法做正規化,只能從 j=1 開始。例如 \(f_𝜃(x) = 𝜃_0 + 𝜃_1x + 𝜃_2x^2\) 中,m 為 2,正規化的參數對象為 \(𝜃_1, 𝜃_2\)。\(𝜃_0\) 被稱為 bias 項


𝜆 是決定對於正規化項的影響為正的常數,要自己決定用什麼值。


正規化的效果


先將目標函數分為兩項


\( C(𝜃)= \frac{1}{2} \sum_{i=1}^{n}( y^{(i)} - f_𝜃(x^{(i)})^2 \)


\(R(𝜃) = \frac{𝜆}{2} \sum_{j=1}^{m} 𝜃_j^2​\)


因為 C(𝜃) 假設是任意一個曲線,R(𝜃) 是二次函數,任意假設它為 \( \frac{1}{2} 𝜃_1^2​\) ,是通過原點的二次函數



正規化後,\(𝜃_1\) 最小值會往原點靠近




\( f_𝜃(x) = 𝜃_0 + 𝜃_1x+ 𝜃_2x^2 \) 是二次曲線,但如果 \( 𝜃_2 \) 為 0,變成一次直線,就簡化了模型


𝜆 的大小,決定正規化的影響,如果 \( 𝜆=0 \) 就等於沒有用到正規化


分類的正規化


分類問題是用對數似然函數


\(\log L(𝜃) = \log \prod _{i=1}^{n} P( y^{(i)} = 1|x^{(i)} )^{y^{(i)}} P( y^{(i)} = 0|x^{(i)} )^{1-y^{(i)}}​\)


正規化就是再加上 R(𝜃) 的部分,另外因為對數似然指數,原本是要最大化,為了轉換為最小化,加上負號


\(\log L(𝜃) = - \log \prod _{i=1}^{n} P( y^{(i)} = 1|x^{(i)} )^{y^{(i)}} P( y^{(i)} = 0|x^{(i)} )^{1-y^{(i)}} + \frac{𝜆}{2} \sum_{j=1}^{m}𝜃_j^2\)


正規化後的微分


因為 \( E(𝜃) = C(𝜃) + R(𝜃) ​\) ,就分別對 C(𝜃), R(𝜃) 進行偏微分


\(\frac{𝜕C(𝜃)}{𝜕𝜃_j} = \sum_{i=1}^{n}( f_𝜃(x^{(i)} )-y^{(i)} )x_j^{(i)} ​\)


\(R(𝜃) = \frac{𝜆}{2} \sum_{j=1}^{m} 𝜃_j^2 = \frac{𝜆}{2}𝜃_1^2 + \frac{𝜆}{2}𝜃_2^2 + \cdots + \frac{𝜆}{2}𝜃_m^2 \)


\( \frac{𝜕R(𝜃)}{𝜕𝜃_j} = 𝜆𝜃_j​\)


因此參數更新式就改為


\(𝜃_j := 𝜃_j - 𝜂 ( \sum_{i=1}^{n}( f_𝜃(x^{(i)} )-y^{(i)} )x_j^{(i)} + 𝜆𝜃_j )\)


關於 \(𝜃_0​\) 的部分,無法處理正規化,因為 R(𝜃) 以 \(𝜃_0​\) 微分後變成 0




邏輯迴歸是類似的過程


E(𝜃) = C(𝜃) + R(𝜃)


\(\log L(𝜃) = - \log \prod _{i=1}^{n} P( y^{(i)} = 1|x^{(i)} )^{y^{(i)}} P( y^{(i)} = 0|x^{(i)} )^{1-y^{(i)}} + \frac{𝜆}{2} \sum_{j=1}^{m}𝜃_j^2\)


\(R(𝜃) = \frac{𝜆}{2} \sum_{j=1}^{m} 𝜃_j^2\)


微分後,因為 R(𝜃) 以 \(𝜃_0​\) 微分後變成 0


\(𝜃_0:= 𝜃_0 - 𝜂 ( \sum_{i=1}^{n}( f_𝜃(x^{(i)}) - y^{(i)} )x_j^{(i)} ) ​\)


\(𝜃_j := 𝜃_j - 𝜂 ( \sum_{i=1}^{n}( f_𝜃(x^{(i)}) - y^{(i)} )x_j^{(i)} + 𝜆𝜃_j ) \quad\quad\quad j > 0 \)




正規化不是只有一種,目前都是 L2 正規化。 L1正規化,用在判斷非必要的參數,會變成0,用來減少變數的數量。L1 正規化是要減少不必要的變數,L2 正規化是要抑制變數的影響。


學習曲線


乏適 Underfitting


跟 Overfitting 相反的狀況是 Underfitting,也就是找不到適合學習資料的模型。


光只有查看精度,沒辦法判斷是 overfitting 或是 underfitting


以這個圖形為例,看起來學習資料是二次曲線,很難找到一條一次函數的直線,適合這些學習資料,隨著學習資料變多,模型的精度會一直下降



使用學習資料數量較少的模型,去預測未知資料會比較困難,精度會比較低。如果學習資料越多,精度會增加。


將學習用資料與測試用資料數量,及精度的對照,繪製圖表。


如出現高偏差的狀況,對學習資料數量增加,但精度降低,測試資料增加,精度會增加,且兩個精度會越來越接近。這種狀況就是 underfitting



如果發現有高方差的狀況,就是對學習資料數量增加,但精度一直維持很高,測試資料增加,卻沒辦法增加精度。這種狀況就是 overfitting



這種將資料個數與精度繪製的圖表,就稱為學習曲線。


References


練好機器學習的基本功:用Python進行基礎數學理論的實作

2019/8/25

機器學習_分類

如果要判斷矩形是橫向或是縱向,通常可以直接看出結果。


可以將所有矩形左下角跟原點對齊,以右上角的座標位置,來判斷橫或縱向。分類就是要找出一條分類的線,將兩種矩形分類。



這條線,是將權重向量 (weight) 視為法線的直線,因為互相垂直,所以內積為 0。


\(w \cdot x = \sum_{i=1}^{n} w_ix_i = 0\)


如果有兩個維度,且 \(weight = (1,1)\)


\(w \cdot x = w_1x_1 + w_2x_2 = 1 \cdot x_1 + 1 \cdot x_2 = x_1 + x_2 = 0\)



內積也可以用另一個式子


\(w \cdot x = |w| \cdot |x| \cdot cos𝜃​\)


因為內積為 0,表示 \(cos𝜃 =0\),也就是 \(𝜃=90^o \) 或 \(𝜃=270^o\)


一般來說是要用機器學習,找出權重向量,得到跟該向量垂直的直線,再透過該直線進行分類。


感知器


感知器是一個可以接受多個輸入,並對每一個值,乘上權重再加總,輸出得到的結果的模型。



  • 準備機器學習的資料

矩形大小 形狀 \(x_1\) \(x_2\) \(y\)
80 x 150 縱向 80 150 -1
60 x 110 縱向 60 110 -1
160 x 50 橫向 160 50 1
125 x 30 橫向 125 30 1

識別函數 \(f_w(x)​\) 就是給定向量 \(x​\) 後,回傳 1 或 -1 的函數,用來判斷橫向或縱向。


\(f_w(x) = \left\{\begin{matrix} 1 \quad (w \cdot x \geq 0) \\ -1 \quad (w \cdot x < 0) \end{matrix}\right.​\)


回到 \(w \cdot x = |w| \cdot |x| \cdot cos𝜃\) 這個式子,如果內積為負數,那麼表示 \( 90^o < 𝜃 < 270^o\)


內積是用來表示向量之間相似程度,如果是正數,就是相似,0 是直角,負數代表不相似。


  • 權重更新式

\(w := \left\{\begin{matrix} w + y^{(i)}x^{(i)} \quad (f_w(x) \neq y^{(i)}) \\ w \quad \quad \quad \quad (f_w(x) = y^{(i)}) \end{matrix}\right.​\)


i 是學習資料的索引,這個權重更新是針對所有學習資料重複處理,用來更新權重向量。上面的部分,意思是藉由識別函數進行分類失敗時,才要去更新權重向量。


權重向量是用隨機的值初始化的



因為一開始識別函數的結果為 -1,而學習資料 \(y^{(1)} =1\),所以要更新權重向量


\( w + y^{(1)}x^{(1)} = w + x^{(1)}​\)



運用向量加法得到新的 w,而新的 w 的法線,會讓 (125, 30) 跟 w 向量在同一側。


線性可分離


感知器有個缺點,只能用來解決線性可分離的問題,以下這樣的問題,無法用一條線去分類。



所以如果要處理圖片分類,就沒辦法以線性的方式處理。


上面例子是單層的感知器,多層感知器,就是神經網路。


另外有一種方法,可用在線性不可分離的問題上:邏輯迴歸 Logistic Regression


邏輯迴歸 Logistic Regression


這種方法是將分類用機率來思考。以一開始矩形橫向或縱向的例子,這邊假設橫向為 1,縱向為 0。


S 型函數


前面的回歸有提到 \(f_𝜃(x) = 𝜃^T x​\) 這個函數,可用最速下降法或隨機梯度下降法學習 𝜃,然後用 𝜃 求得未知資料 x 的輸出值。


這邊需要的函數為,其中 \(exp(-𝜃^Tx) = e^{-𝜃^T x}\) , e 為自然對數的底數 Euler's number (2.71828)


\(f_𝜃(x) = \frac{1}{1 + exp(-𝜃^Tx)}​\)


會稱為 S 函數是因為如果將 \(𝜃^Tx\) 設為橫軸,\(f_𝜃(x)\) 設定為縱軸,會出現這樣的圖形



S函數的特徵是 當 \(𝜃^Tx = 0​\) 會得到 \(f_𝜃(x) = 0.5​\),且 \(0< f_𝜃(x) <1 ​\)


當作機率處理的原因是 \(0< f_𝜃(x) <1 \)


決策邊界


將未知資料 x 屬於橫向的機率設為 \(f_𝜃(x)\) ,用條件機率的方式描述為


\(P( y = 1 | x) = f_𝜃(x)​\)


當給予資料 x 時,y=1的機率為 \(f_𝜃(x)\) 。如果計算後的結果,機率為 0.7,就表示矩形為橫向的機率為 0.7。以 0.5 為閥值,判斷是不是橫向。


\(y = \left\{\begin{matrix} 1 \quad (f_𝜃(x) \geq 0.5) \\ 0 \quad (f_𝜃(x) < 0.5) \end{matrix}\right.​\)


回頭看 S 函數,當 \(f_𝜃(x) \geq 0.5​\),也就是 \(𝜃^T x >0​\) ,就判定為橫向。可將判斷式改為


\(y = \left\{\begin{matrix} 1 \quad (𝜃^T x \geq 0) \\ 0 \quad (𝜃^T x < 0) \end{matrix}\right.\)




任意選擇一個 𝜃 為例子,\(x_1\) 是橫長, \(x_2\) 是高


\(𝜃= \left[ \begin{matrix} 𝜃_0 \\ 𝜃_1 \\𝜃_2 \end{matrix} \right] = \left[ \begin{matrix} -100 \\ 2 \\ 1 \end{matrix} \right] , x= \left[ \begin{matrix} 1 \\ x_1 \\x_2 \end{matrix} \right] ​\)


\(𝜃^T x = -100 \cdot 1 +2x_1 +x_2 \geq 0 \)


\( x_2 \geq -2x_1 +100 \) 就表示分類為橫向


以圖形表示



以 \(𝜃^Tx =0\) 這條直線為邊界線,就能區分橫向或縱向,這條線就是決策邊界。但實際上,這個任意選擇的 𝜃 並不能正確的進行分類,因此為了求得正確的 𝜃 ,就要定義目標函數,進行微分,以求得正確的參數 𝜃 ,這個方法就稱為邏輯迴歸。


似然函數


現在要找到 𝜃 的更新式


一開始將 x 為橫向的機率 \(P( y = 1 | x) ​\) 定義為 \(f_𝜃(x)​\) ,根據這個定義,學習資料 \(y​\) 跟 \(f_𝜃(x)​\) 的關係,最佳的狀況是 \(y=1, f_𝜃(x)=1​\) , \(y=0, f_𝜃(x)=0​\) ,但還要改寫為


  • \( y=1 ​\) 時,要讓機率 \(P( y = 1 | x) ​\) 是最大,判定為橫向
  • \( y=0 ​\) 時,要讓機率 \(P( y = 0| x) ​\) 是最大,判定為縱向

矩形大小 形狀 \(y\) 機率
80 x 150 縱向 0 要讓機率 $P( y = 0
60 x 110 縱向 0 要讓機率 $P( y = 0
160 x 50 橫向 1 要讓機率 $P( y = 1
125 x 30 橫向 1 要讓機率 $P( y = 1

因為所有學習資料互相獨立沒有關聯,整體的機率就是全部的機率相乘


\(L(𝜃) = P( y^{(1)} = 0|x^{(1)} ) P( y^{(2)} = 0|x^{(2)} ) P( y^{(3)} = 1|x^{(3)} ) P( y^{(4)} = 1|x^{(4)} ) \)


這個式子可改寫為


\(L(𝜃) = \prod _{i=1}^{n} P( y^{(i)} = 1|x^{(i)} )^{y^{(i)}} P( y^{(i)} = 0|x^{(i)} )^{1-y^{(i)}}​\)




如果假設 \(y^{(i)} = 1​\)


\(P( y^{(i)} = 1|x^{(i)} )^1 P( y^{(i)} = 0|x^{(i)} )^0 = P( y^{(i)} = 1|x^{(i)} )​\)


如果假設 \(y^{(i)} =0 \)


\(P( y^{(i)} = 1|x^{(i)} )^0 P( y^{(i)} = 0|x^{(i)} )^1 = P( y^{(i)} = 0|x^{(i)} )​\)




目標函數 \(L(𝜃)\) 就稱為似然函數 Likelihood,就是要找到讓 \(L(𝜃)\) 最大的參數 𝜃


對數似然函數


因為機率都小於 1,機率的乘積會不斷變小,在程式設計會產生精確度的問題。所以加上 log


\(\log L(𝜃) = \log \prod _{i=1}^{n} P( y^{(i)} = 1|x^{(i)} )^{y^{(i)}} P( y^{(i)} = 0|x^{(i)} )^{1-y^{(i)}}​\)


因為 log 是單調遞增函數,因此不會影響到結果。換句話說,要讓 \(L(𝜃)​\) 最大化,跟要讓 \(logL(𝜃)​\) 最大化是一樣的。


\(\begin{equation}
\begin{split}
\log L(𝜃) &= \log \prod _{i=1}^{n} P( y^{(i)} = 1|x^{(i)} )^{y^{(i)}} P( y^{(i)} = 0|x^{(i)} )^{1-y^{(i)}}\\
&=\sum_{i=1}^{n}( \log P( y^{(i)} = 1|x^{(i)} )^{y^{(i)}} + log P( y^{(i)} = 0|x^{(i)} )^{1-y^{(i)}} ) \\
&=\sum_{i=1}^{n}( {y^{(i)}} \log P( y^{(i)} = 1|x^{(i)} ) + ({1-y^{(i)}}) \log P( y^{(i)} = 0|x^{(i)} ) ) \\
&=\sum_{i=1}^{n}( {y^{(i)}} \log P( y^{(i)} = 1|x^{(i)} ) + ({1-y^{(i)}}) \log (1-P( y^{(i)} = 1|x^{(i)} )) ) \\
&=\sum_{i=1}^{n}( {y^{(i)}} \log f_𝜃( x^{(i)} ) + ({1-y^{(i)}}) \log (1- f_𝜃(x^{(i)} )) ) \\
\end{split}
\end{equation}​\)


  • \(\log(ab) = \log a + \log b\)
  • \(\log a^b = b \log a​\)
  • 因為只考慮 \(y=1\) 或 \(y=0\),所以 \(P( y^{(i)} = 0|x^{(i)} ) = 1 - P( y^{(i)} = 1|x^{(i)} )\)

似然函數的微分


邏輯迴歸,就是將這個對數似然函數當作目標函數使用


\( \log L(𝜃) =\sum_{i=1}^{n}( {y^{(i)}} \log f_𝜃( x^{(i)} ) + ({1-y^{(i)}}) \log (1- f_𝜃(x^{(i)} )) ) ​\)


要將這個函數,個別針對參數 \(𝜃_j\) 進行偏微分


同樣利用合成函數的微分方法


\( u = \log L(𝜃)​\)


\(v = f_𝜃 (x) = \frac{1}{1 + exp(-𝜃^Tx)}​\)


然後


\(\frac{𝜕E}{𝜕𝜃_j} = \frac{𝜕u}{𝜕v} \frac{𝜕v}{𝜕𝜃_j}​\)




先計算第一項


因為 \(\log (v)\) 的微分是 \(\frac{1}{v}\)


而 \( \log (1-v)\) 的微分為


\( s =1-v ​\)


\( t = \log (s) \)


\( \frac{dt}{dv} = \frac{dt}{ds} \cdot \frac{ds}{dv} = \frac{1}{s} \cdot -1 = - \frac{1}{1-v} \)


所以


\( \frac{𝜕u}{𝜕v} = \frac{𝜕}{𝜕v} \sum_{i=1}^{n}( {y^{(i)}} \log (v) + ({1-y^{(i)}}) \log (1- v ) ) = \sum_{i=1}^{n} ( \frac{y^{(i)}}{v} - \frac{1- y^{(ii)} }{1-v} )\)




然後將 \(v\) 以 \(𝜃_j\) 微分


\(\frac{𝜕v}{𝜕𝜃_j} = \frac{𝜕}{𝜕𝜃_j} \frac{1}{ 1+ exp(-𝜃^Tx)} ​\)


因為 \(f_𝜃 (x)\) 是 S 型函數,且已知 S 型函數的微分為


\( \frac{d𝜎(x)}{dx} = 𝜎(x) (1-𝜎(x)) \)


利用合成函數的微分方法


\( z = 𝜃^T x\)


\(v = f_𝜃 (x) = \frac{1}{1 + exp(-z)}\)


\(\frac{𝜕v}{𝜕𝜃_j} = \frac{𝜕v}{𝜕z} \frac{𝜕z}{𝜕𝜃_j} ​\)


前面的部分


\( \frac{𝜕v}{𝜕z} = v(1-v) ​\)


後面的部分


\( \frac{𝜕z}{𝜕𝜃_j} = \frac{𝜕}{𝜕𝜃_j} 𝜃^Tx = \frac{𝜕}{𝜕𝜃_j} (𝜃_0x_0 +𝜃_1x_1 +\cdots + 𝜃_nx_n ) = x_j ​\)


所以


\(\frac{𝜕v}{𝜕𝜃_j} = \frac{𝜕v}{𝜕z} \frac{𝜕z}{𝜕𝜃_j} = v(1-v) x_j ​\)




\(\begin{equation}
\begin{split}
\frac{𝜕u}{𝜕𝜃_j} &= \frac{𝜕u}{𝜕v} \frac{𝜕v}{𝜕𝜃_j} \\
& = \sum_{i=1}^{n}( \frac{y^{(i)}}{v} - \frac {1 - y^{(i)}}{1-v} ) \cdot v(1-v) \cdot x_j^{(i)} \\
& = \sum_{i=1}^{n}( y^{(i)}(1-v) - (1-y^{(i)})v )x_j^{(i)} \\
& = \sum_{i=1}^{n}( y^{(i)} -v )x_j^{(i)} \\
& = \sum_{i=1}^{n}( y^{(i)} - f_𝜃(x^{(i)}) )x_j^{(i)} \\
\end{split}
\end{equation}​\)


先前最小化,是要往微分後的結果的正負符號相反方向移動。但現在要最大化,所以要往微分後的結果的正負符號相同方向移動。


\( 𝜃_j := 𝜃_j + 𝜂 \sum_{i=1}^{n}( y^{(i)} - f_𝜃(x^{(i)}) )x_j^{(i)} ​\)


也能配合多元迴歸,改寫成這樣


\(𝜃_j := 𝜃_j - 𝜂 \sum_{i=1}^{n}( f_𝜃(x^{(i)}) - y^{(i)} )x_j^{(i)}\)


線性不可分離


線性不可分離的問題不能直線,但可嘗試用曲線。


例如,將 \(x_1^2\) 加入學習資料


\(𝜃= \left[ \begin{matrix} 𝜃_0 \\ 𝜃_1 \\𝜃_2 \\𝜃_3 \end{matrix} \right] , x= \left[ \begin{matrix} 1 \\ x_1 \\x_2 \\x_1^2\end{matrix} \right] ​\)


然後


\( 𝜃^Tx = 𝜃_0+𝜃_1x_1+𝜃_2x_2+𝜃_3x_1^2 ​\)


假設


\(𝜃= \left[ \begin{matrix} 𝜃_0 \\ 𝜃_1 \\𝜃_2 \\𝜃_3 \end{matrix} \right] = \left[ \begin{matrix} 0 \\ 0 \\1 \\-1 \end{matrix} \right] \)


因為 \(𝜃^Tx \geq 0​\)


\( 𝜃^Tx = 𝜃_0+𝜃_1x_1+𝜃_2x_2+𝜃_3x_1^2 = x_2 - x_1^2 \geq 0 ​\)


得到方程式 \( x_2 \geq x_1^2 ​\)



現在的決策邊界變成曲線,因為參數 𝜃 是任意選擇的,所以資料沒有被正確地分類。


可以增加次方數,得到複雜的形狀的決策邊界。


另外還有 SVM (支援向量機) 的分類演算法,多元分類處理方法。


References


練好機器學習的基本功:用Python進行基礎數學理論的實作

2019/8/18

機器學習_線性迴歸


機器學習擅長


  1. 回歸 regression


    將連續性的資料進行觀察,用以預測未來的結果。例如股價、身高、體重

  2. 分類 classification


    收集既有的資料,進行訓練,根據訓練結果預測新資料的分類。例如:垃圾郵件判斷、手寫數字辨識

  3. 分群 clustering


    根據資料進行分群,但跟分類不同,分類的訓練資料已經有標記結果,要用來分群的資料並沒有群組的標記。例如:根據學測成績,進行文理組分群


使用具有標記的資料進行機器學習,稱為監督式學習。


使用不具有標記的資料進行機器學習,稱為非監督式學習。


線性迴歸(Linear regression)


當我們取得原始量測資料時,如果在平面座標上標記這些量測點,會感覺到這些點之間,可以畫出最接近這些點的一條直線方程式,線性回歸方法,可以找到這樣的方程式,未來就可以根據這個方程式,預測數值。



從圖形看起來,我們可找到一條「最接近」所有紅色觀測值的直線,以直線方程式 \({f(x)=ax+b}\) 表示這條直線,我們要做的就是找到一個方法,確定 a 與 b 的值,未來就可以利用這個方程式預測數據。在統計中,為了因應未來可能有很多未知數的問題,改以這樣的寫法:


\({f(x)=𝜃_0+𝜃_1x}\)


最小平方法


假設目前有這些數據,當我們任意找 \({𝜃_0 =1, 𝜃_1=2}​\) 時,f(x) 跟實際上的 y 之間有誤差。所以要找到適當的參數,讓 f(x) 與 y 之間的誤差最小,當然如果誤差為 0 是最好的。


x y f(x)
58 374 117
70 385 141
81 375 163
84 401 169

定義誤差函數為


\({E(𝜃)= \frac{1}{2} \sum_{i=1}^{n}( y^{(i)} - f_𝜃(x^{(i)})^2 }\)


  • \(x^{(i)}, y^{(i)}\) 分別是第 i 項的 x 與 y,例如 \(x^{(1)} = 58, y^{(1)}=374\)
  • \({(y^{(i)} - f_𝜃(x^{(i)}) }\) 是誤差值,但因為誤差有可能是負數,所以就用平方,轉成正數
  • 將所有誤差的平方加總後,為了微分計算方便,就在前面再乘上 1/2。因為乘上任意的正數,只會讓圖形橫向壓扁,但不會改變最小值的位置。任意的正數都可以,選擇 1/2 是因為後面的例子,f(x) 是二次函數的關係。

在讓 E(𝜃) 最小的狀況下,找到的 \(𝜃_0, 𝜃_1​\) 就是最小平方法


最速下降法


剛剛要讓 E(𝜃) 最小的狀況下,必須不斷地找到不同的 \(𝜃_0, 𝜃_1​\),這個計算很麻煩,可用微分來解決,因為微分就是在找函數切線斜率的變化。


例如 \(f(x) = (x-1)^2\) ,微分後 \(f'(x) = 2x-2\)


x f'(x) 的正負 f(x) 遞增或遞減
\(x < 1\) \(-\) 遞減
\(x=0\) 0
\(x>1\) \(+\) 遞增

f(x) 的圖形如下,當 x 由 3 往 1 逼近,f(x) 就越來越小,另外當 x 由 -1 往 1 逼近,f(x) 也會越來越小



意思就是說,只要 x 往導函數 (微分) 的反方向移動,就函數值會往最小值移動。


最速下降法(梯度下降法) 就是定義為


\(x := x - 𝜂 \frac{d}{dx}f(x)​\)


以實際的數字為例,當 \(𝜂 = 1, x = 3​\),x 會在 3 與 -1 之間往返


\(x := 3-1(2*3-2) = -1​\)


\(x := -1-1(2*(-1)-2) = 3​\)


當 \(𝜂 = 0.1, x = 3\),x 會往最小值逼近


\(x := 3-0.1(2*3-2) = 2.6​\)


\(x := 2.6-0.1(2*2.6-2) = 2.3​\)


\(x := 2.3-0.1(2*2.3-2) = 2.1​\)


當 𝜂 越大,x 就會往返,當 𝜂 越小,x 會往最小值逼近




回到剛剛的 誤差函數


\({E(𝜃)= \frac{1}{2} \sum_{i=1}^{n}(y^{(i)} - f_𝜃(x^{(i)})^2 }\)


因為 \(f_𝜃(x^{(i)})\) 是 \({f(x)=𝜃_0+𝜃_1x}\) ,有兩個未知的參數 \(𝜃_0, 𝜃_1\) ,要改用偏微分找最小值。


\(𝜃_0 := 𝜃_0 - 𝜂 \frac{𝜕E}{𝜕𝜃_0}​\)


\(𝜃_1 := 𝜃_1 - 𝜂 \frac{𝜕E}{𝜕𝜃_1}​\)


因 E(𝜃) 裡面有 \(f_𝜃(x)​\),而 \(f_𝜃(x)​\) 裡面有 𝜃


\(u = E(𝜃)​\)


\(v = f_𝜃(x)​\)


然後用合成函數的方式,計算微分


\(\frac{𝜕E}{𝜕𝜃_0} = \frac{𝜕u}{𝜕v} \frac{𝜕v}{𝜕𝜃_0}​\)


其中,前面的部分


\( \frac{𝜕u}{𝜕v} = \frac{𝜕}{𝜕v}( \frac{1}{2} \sum_{i=1}^{n}(y^{(i)} - v)^2 ) = \frac{1}{2} \sum_{i=1}^{n}\frac{𝜕}{𝜕v}(y^{(i)} - v)^2 = \frac{1}{2}\sum_{i=1}^{n}( -2y^{(i)} +2v ) = \sum_{i=1}^{n}( v-y^{(i)} )\)


後面的部分


\(\frac{𝜕v}{𝜕𝜃_0} = \frac{𝜕}{𝜕𝜃_0}( 𝜃_0 + 𝜃_1x ) = 1​\)


所以


\(\frac{𝜕E}{𝜕𝜃_0} = \sum_{i=1}^{n}( f_𝜃(x^{(i)})-y^{(i)} )\)


另外對 \(𝜃_1​\) 微分,可得到


\(\frac{𝜕v}{𝜕𝜃_1} = \frac{𝜕}{𝜕𝜃_1}( 𝜃_0 + 𝜃_1x ) = x\)


\(\frac{𝜕E}{𝜕𝜃_1} = \sum_{i=1}^{n}( f_𝜃(x^{(i)} )-y^{(i)} )x^{(i)} \)




最後


\(𝜃_0 := 𝜃_0 - 𝜂 \sum_{i=1}^{n}( f_𝜃(x^{(i)} )-y^{(i)} )\)


\(𝜃_1 := 𝜃_1 - 𝜂 \sum_{i=1}^{n}( f_𝜃(x^{(i)} )-y^{(i)} )x^{(i)}​\)


用這個方法,就可以找出正確的 \(𝜃_0, 𝜃_1​\)


多項式回歸


一開始,我們假設數據的模型是線性的,所以使用一次函數,但也可能用二次或更高次的函數來定義 \(f_𝜃(x)​\),會更貼近原本的數據模型


\(f_𝜃(x) = 𝜃_0 + 𝜃_1x + 𝜃_2x^2​\)


\(f_𝜃(x) = 𝜃_0 + 𝜃_1x + 𝜃_2x^2 + \dots +𝜃_nx^n​\)


回到剛剛的問題,要對 \(𝜃_2​\) 進行偏微分


對 \(𝜃_1​\) 微分,可得到


\(\frac{𝜕v}{𝜕𝜃_2} = \frac{𝜕}{𝜕𝜃_2}( 𝜃_0 + 𝜃_1x +𝜃_2x^2 ) = x^2\)


\(\frac{𝜕E}{𝜕𝜃_2} = \sum_{i=1}^{n}( f_𝜃(x^{(i)} )- y^{(i)} )(x^{(i)} )^2\)


多元回歸


目前解決的問題,都只有一個變數 x,但大多數的問題,都是有兩個以上的變數。例如廣告的點擊率,可能會受廣告費、顯示位置、顯示大小 等原因影響。


\(f_𝜃(x_1, x_2, x_3)=𝜃_0+𝜃_1x_1+𝜃_2x_2+𝜃_3x_3\)


當變數有 n 個,可改用向量的方式表示


\(𝜃= \left[ \begin{matrix} 𝜃_0 \\ 𝜃_1 \\𝜃_2 \\. \\. \\𝜃_n \end{matrix} \right] x= \left[ \begin{matrix} x_1 \\ x_2 \\x_3 \\. \\. \\x_n \end{matrix} \right] ​\)


但因為 𝜃 跟 x 個數不同,不容易計算,就再加上一項 \(x_0 =1​\)


\(𝜃= \left[ \begin{matrix} 𝜃_0 \\ 𝜃_1 \\𝜃_2 \\. \\. \\𝜃_n \end{matrix} \right] x= \left[ \begin{matrix} x_0 \\ x_1 \\x_2 \\. \\. \\x_n \end{matrix} \right] \)


將 𝜃 變成轉置矩陣後,再跟 x 相乘,就會是剛剛的 \(f_𝜃(x)\)


\(𝜃^Tx = 𝜃_0x_0+𝜃_1x_1+ \dots + 𝜃_nx_n = f_𝜃(x) \)


變成向量後,再用剛剛合成函數偏微分的方法


\(u = E(𝜃)​\)


\(v = f_𝜃(x)​\)


\(\frac{𝜕u}{𝜕𝜃_j} =\frac{𝜕E}{𝜕𝜃_j} = \frac{𝜕u}{𝜕v} \frac{𝜕v}{𝜕𝜃_j}\)


前面的部分一樣,後面的部分


\(\frac{𝜕v}{𝜕𝜃_j} = \frac{𝜕}{𝜕𝜃_j}( 𝜃^Tx ) = \frac{𝜕}{𝜕𝜃_j}( 𝜃_0x_0+𝜃_1x_1+\dots+𝜃_nx_n )= x_j\)


第 j 項參數的定義為


\(𝜃_j := 𝜃_j - 𝜂 \sum_{i=1}^{n}( f_𝜃(x^{(i)} )-y^{(i)} )x_j^{(i)}\)


當變數增加,計算量變大,用最速下降法會導致計算速度變慢,可用隨機梯度下降法改進。


最速下降法除了有計算速度慢的問題,還有可能陷入局部解的問題,像以下的函數圖形中,不同的起點,可能會找到局部最小值。



隨機梯度下降法


在多元迴歸中,第 j 項參數的定義為


\(𝜃_j := 𝜃_j - 𝜂 \sum_{i=1}^{n}( f_𝜃(x^{(i)} )-y^{(i)} )x_j^{(i)}​\)


但因為用到所有的資料的誤差,計算量太大,隨機梯度下降法式隨機選擇一項學習資料,套用在參數的更新上,例如選擇第 k 項。


\(𝜃_j := 𝜃_j - 𝜂 ( f_𝜃(x^{(k)} )-y^{(k)} )x_j^{(k)}​\)


原本最速下降法用來更新一次參數的時間,隨機梯度下降法可更新 n 次參數。因為是隨機選擇學習資料,不會陷入局部解的問題。




另外也有隨機選擇 m 筆學習資料的方法,也稱為小量批次資料法,假設 m 筆資料的集合為 K


\(𝜃_j := 𝜃_j - 𝜂 \sum_{k𝜖K} ( f_𝜃(x^{(k)} )-y^{(k)} )x_j^{(k)}​\)


References


練好機器學習的基本功:用Python進行基礎數學理論的實作

2019/8/12

erlang lager with date in log filename

erlang lager 預設是以設定中的 filename 加上 .1 .2 的 postfix 作為 logfile rotate 的依據,但通常在使用 logfile,會希望直接在 logfile 看到產生該 log 的日期,這時需要使用 Custom Log Rotation 的功能,自己撰寫 log_rotator。


首先我們先找到 lager 原始程式碼中預設的 lagerrotatordefault.erl,先複製成 mylagerlog_rotator,然後修改裡面的程式碼。


-module(my_lager_log_rotator).

-include_lib("kernel/include/file.hrl").

-behaviour(lager_rotator_behaviour).

-export([
  create_logfile/2, open_logfile/2, ensure_logfile/4, rotate_logfile/2
]).

create_logfile(Name, Buffer) ->
  {{Y, M, D}, {H, _, _}} = calendar:now_to_local_time(os:timestamp()),
  DateHour =  {Y, M, D, H},
  FileName = filename(Name, DateHour, 1),
  file:delete(Name),
  file:make_symlink(filename:absname(FileName), Name),
  open_logfile(Name, Buffer).

open_logfile(Name, Buffer) ->
  case filelib:ensure_dir(Name) of
    ok ->
      Options = [append, raw] ++
        case  Buffer of
          {Size, Interval} when is_integer(Interval), Interval >= 0, is_integer(Size), Size >= 0 ->
            [{delayed_write, Size, Interval}];
          _ -> []
        end,
      case file:open(Name, Options) of
        {ok, FD} ->
          case file:read_file_info(Name) of
            {ok, FInfo} ->
              Inode = FInfo#file_info.inode,
              {ok, {FD, Inode, FInfo#file_info.size}};
            X -> X
          end;
        Y -> Y
      end;
    Z -> Z
  end.


ensure_logfile(Name, FD, Inode, Buffer) ->
  case file:read_link(Name) of
    {ok, _} ->
      lager_ensure_logfile(Name, FD, Inode, Buffer);
    _ ->
      create_logfile(Name, Buffer)
  end.


lager_ensure_logfile(Name, undefined, _Inode, Buffer) ->
  open_logfile(Name, Buffer);
lager_ensure_logfile(Name, FD, Inode, Buffer) ->
  case file:read_file_info(Name) of
    {ok, FInfo} ->
      Inode2 = FInfo#file_info.inode,
      case Inode == Inode2 of
        true ->
          {ok, {FD, Inode, FInfo#file_info.size}};
        false ->
          %% delayed write can cause file:close not to do a close
          _ = file:close(FD),
          _ = file:close(FD),
          case open_logfile(Name, Buffer) of
            {ok, {FD2, Inode3, Size}} ->
              %% inode changed, file was probably moved and
              %% recreated
              {ok, {FD2, Inode3, Size}};
            Error ->
              Error
          end
      end;
    _ ->
      %% delayed write can cause file:close not to do a close
      _ = file:close(FD),
      _ = file:close(FD),
      case open_logfile(Name, Buffer) of
        {ok, {FD2, Inode3, Size}} ->
          %% file was removed
          {ok, {FD2, Inode3, Size}};
        Error ->
          Error
      end
  end.
%%
%%%% renames failing are OK
%%rotate_logfile(File, 0) ->
%%  %% open the file in write-only mode to truncate/create it
%%  case file:open(File, [write]) of
%%    {ok, FD} ->
%%      file:close(FD),
%%      ok;
%%    Error ->
%%      Error
%%  end;
%%rotate_logfile(File0, 1) ->
%%  File1 = File0 ++ ".0",
%%  _ = file:rename(File0, File1),
%%  rotate_logfile(File0, 0);
%%rotate_logfile(File0, Count) ->
%%  File1 = File0 ++ "." ++ integer_to_list(Count - 2),
%%  File2 = File0 ++ "." ++ integer_to_list(Count - 1),
%%  _ = file:rename(File1, File2),
%%  rotate_logfile(File0, Count - 1).
%%

rotate_logfile(Name, _Count) ->
  case file:read_link(Name) of
    {ok, LinkedName} ->
      case filelib:file_size(LinkedName) of
        0 ->
          %% if the files size is zero, it is removed
          catch file:delete(LinkedName);
        _ ->
          void
      end;
    _ ->
      void
  end,
  {ok, {FD, _, _}} = create_logfile(Name, []),
  file:close(FD).

%% @doc Create name of a new file
%% @private
filename(BaseFileName, DateHour, Branch) ->
  FileName = lists:append([BaseFileName,
    suffix(DateHour, false), ".", integer_to_list(Branch)
  ]),
  case filelib:is_file(FileName) of
    true ->
      filename(BaseFileName, DateHour, Branch + 1);
    _ ->
      FileName
  end.

%% @doc Zero-padding number
%% @private
zeropad(Num, MinLength) ->
  NumStr = integer_to_list(Num),
  zeropad_str(NumStr, MinLength - length(NumStr)).
zeropad_str(NumStr, Zeros) when Zeros > 0 ->
  zeropad_str([$0 | NumStr], Zeros - 1);
zeropad_str(NumStr, _) ->
  NumStr.

%% @doc Create a suffix
%% @private
suffix({Y, M, D, H}, WithHour) ->
  YS = zeropad(Y, 4),
  MS = zeropad(M, 2),
  DS = zeropad(D, 2),
  HS = zeropad(H, 2),
  case WithHour of
    true ->
      lists:flatten([$., YS, MS, DS, $., HS]);
    _ ->
      lists:flatten([$., YS, MS, DS])
  end.

將 mylagerlogrotator 套用在 lager 的設定檔的 lagerfile_backend 中,{rotator, my_lager_log_rotator}


[
  {lager, [
    {log_root, "./log"},
    {crash_log, "crash.log"},
    {error_logger_redirect, false},
    {colored, true},
    {colors, [
      {debug,     "\e[0;36m" },
      {info,      "\e[1;37m" },
      {notice,    "\e[1;36m" },
      {warning,   "\e[1;33m" },
      {error,     "\e[1;31m" },
      {critical,  "\e[1;35m" },
      {alert,     "\e[1;44m" },
      {emergency, "\e[1;41m" }
    ]},
    {handlers, [
      {lager_console_backend, [{level, debug}, {formatter, lager_default_formatter},
        {formatter_config, [date, " ", time, color, " ", pid, " ", module, ":", line, " [", severity, "] ", message, "\e[0m\n"]}]},
      {lager_file_backend, [{file, "debug.log"}, {level, debug}, {size, 3000}, {date, "$H00"}, {count, 2},
        {formatter_config, [date, " ", time, " ", pid, " ", module, ":", line, " [", severity, "] ", message, "\n"]}, {rotator, my_lager_log_rotator}]}
    ]}
  ]}
].

現在就會產生像這樣的 logfile


debug.log (symbolic link to debug.log.20190426.2)
debug.log.20190426.1
debug.log.20190426.2

目前還會需要調整的是設定檔中的 count,如果在 logfile 加上日期,count 應該是要代表保留幾天的資料,但目前還是依照 lager 原本的定義,為保留幾個同樣 prefix 檔名的 logfile。


References


erlang lager


leologgerrotator.erl