Note20200605

這是我去雄女研習的時候做的筆記
有點簡陋啦 不過看得懂就沒差
原本是塞在 Hackmd 的我就把它貼過來吧

# 前言

還記得筆者國中做科展的時候
都是用 Excel 裡面的圖表功能
把數據整理出來推測出合理的函數
當時都不知道是如何做到的
這次就要來實作

# 前提

• 使用 Py 及其下套件實作 以讀者具有 Py 語法先備知識為前提寫作
• 需有基本線性規劃及微積分數學基礎
• 你非常清楚自己在寫什麼 code

# 線性回歸

說到 Excel 圖表裡的回歸功能
在數學上叫做線性回歸
當我們給定大量測資時
我們不能保證所有測資可以形成一個漂亮的函數圖形
因此我們就會需要這個功能來做輔助
實際上就是給定大量數據讓 AI 去推演出較為合理的函數圖形
要怎麼做呢 讓我們往下看

# 預測模型

先看圖

圖中的小點是我們給予電腦的資料
也就是有著一堆( x_i , y_i )的群組
而斜直線則是 AI 預測的函數 (L : Y=β₀+β₁x)
可想而知
我們要獲得最符合資料統整結果的直線 L
就需要讓每個測資對應函數的距離達到最小
亦即ε_i=y_i-Y_i須達到最小化

# 找最小值

那我們要如何才能找到最適合的回歸函數呢
換句話說 我們要找到最適合的 β₀和 β₁
其實不難
只要把對於所有測資的ε_i找出來相加再找出最小值就好
首先要處理的第一件事
要先處理掉 ε_i的正負差
這裡要注意一點ε_i在算的時候只要y_i<Y_i
ε_i就會是負的
為了避免掉正負相加會抵消
這裡要搬出一個原則叫最小平方法
先平方再加能有效避免掉正負相間的問題

好的 這時候我們可以得出一個式子

$$L( β_0 , β_1 )=\frac {1} {2}Σ[y_i-(β_0+β_1x)]^2$$

接下來為了讓目標函數具有最小值:
目標函數分別對 β0、β1 取偏微分的值必須為零

• $$\frac {∂L( β_0 , β_1 )} {∂β_0}=Σ[y_i-(β_0+β_1x)]=0$$

• $$\frac {∂L( β_0 , β_1 )} {∂β_1}=Σ[y_i-(β_0+β_1x)]x_i=0$$

大概長這樣 (偷偷幫你省略掉中間一大串化簡步驟)

# 梯度下降

這時候你可能會想
欸？難不成我們直接用公式的方式得到最小值嗎
但我們要如何找到那個所謂的最小值呢
這時候要搬出一個方法 叫梯度下降

地球人都知道
求導函數即為求該點的斜率
而當導函數值為 0 時會有極值
此時我們可以有個大膽的想法
先看看圖

發現了嗎 你得到的斜率會告訴你該往哪走才會到極值 (也就是所謂的谷底)

• 若在函數圖形上切線斜率大於 0 的位置
  則極值在 X 向左移的方向
• 若在函數圖形上切線斜率小於的位置
  則極值在 X 向右移的方向

我們可以讓 x 慢慢移動使其接近谷底
修正的方法為

$$x=x-η\frac{dy}{dx}$$

η 我們稱之為學習率 意即每次修正的幅度
這裡要注意的是

• 學習率 η 若太大可能發生無法收斂
• 學習率 η 若太小可能發生收斂速度太慢

因此我們要讓 AI 去學習 逐步去修正 β₀跟 β₁
可以列出式子

• $$β_0:β_0-η\frac {∂L( β_0 , β_1 )} {∂β_0}=β_0-ηΣ[y_i-(β_0+β_1x)]$$

• $$β_1:β_1-η\frac {∂L( β_0 , β_1 )} {∂β_1}=β_1-ηΣ[y_i-(β_0+β_1x)]x_i$$

此時大致上理論就成形了
讓我們來試做看看

# 實作

理論的部分先緩一下
讓我們把上面那坨東東先打成程式

首先先 define 出我們的目標函數 好讓我們之後驗證可以用

def f(x):
    return b0 + b1 * x

當然 b0 跟 b1 不能白有
這裡我們選擇讓它們隨機開始
(由於後面還會用到 於是這裡是直接使用 numpy 裡面的 rand)

import numpy as np

b0 = np.random.rand()
b1 = np.random.rand()

再來就是要給予我們的資料值
這裡使用矩陣來存 方便我們運算可以一口氣運算

data_x=np.array([裡面是一坨資料])
data_y=np.array([裡面是一坨資料])

接著我們分割資料
讓一部份的資料可以拿來學習
另一部分的資料可以拿來驗算
當然你也可以選擇全部拿來學習 再另行人工驗證也可以
這裡採用的是前者 分成兩半進行

t_x=data_x[:len(data_x)//2]
t_y=data_y[:len(data_x)//2]

v_x=data_x[len(data_x)//2:]
v_y=data_y[len(data_x)//2:]

接著就是訓練了
我們要讓機器透過學習來修正自己的參數
我們姑且讓機器學習個 20000 次
並將學習率設定為 0.0001

ETA=0.0001
t_ct=20000

def train():
    global b0
    global b1
    for i in range(t_ct):
        b0 -= ETA * np.sum((f(t_x) - t_y))
        b1 -= ETA * np.sum((f(t_x) - t_y) * t_x)

Do Re Mi So~
這樣就 AI 就學習了 20000 次了
理論上 β₀跟 β₁就已經被修正到接近目標函數了
接下來就讓我們看看結果吧
我們使用 matplotlib 的圖表功能來顯示我們的結果

import matplotlib.pyplot as plt

def check():
    x = np.linspace(0,30, 2000)
    plt.figure(figsize=(8,6))
    plt.title("train-data",fontsize=20,color='white')
    plt.xticks(range(0,35,5),color='white')
    plt.yticks(range(300,700,100),color='white')
    plt.xlabel("x",fontsize=20,color='white')
    plt.ylabel("y",fontsize=20,color='white')
    plt.plot(x, f(x),label="mod func")
    plt.plot(t_x, t_y, 'go',label="train data")
    plt.plot(v_x, v_y, 'ro',label="valid data")
    plt.legend()
    plt.show()

結果長這樣
漂亮吧

# 優化

這時候出現了個問題
當我們某個特徵範圍過大 而另一個太小
當我們在做梯度下降時
收斂的複雜度會比特徵平均的資料還要來的複雜許多
這時我們可以透過一些簡單的優化來提升效率

# 特徵縮放

如上所述
特徵縮放是將特徵資料按比例縮放
讓資料落在某一特定的區間
去除數據的單位限制
將其轉化為純數值
以便不同單位的特徵能夠進行比較
除了可以優化梯度下降法
還能提高精確度

講起來很複雜 大致上就是高一統計學的兩種方法

1. 標準化
   將所有特徵數據縮放成平均為 0、標準差為單位
2. 歸一化
   將特徵數據按比例縮放到 0 到 1 的區間

做起來一點也不難

# 標準化

$$x'=\frac{data-min(data)}{fullrange(data)}$$

# 歸一化

$$x'=\frac{data-arg(data)}{std(data)}$$

簡單吧 做成程式也不難
只要把 data_x 修改為一下就好
在程式中加入

def fc(x,t):
    return (x-np.min(x))/(np.max(x)-np.min(x)) if t==1 else (x-np.mean(x))/np.std(x)

data_x=fc(data_x)

為了能方便選擇要標準化還是歸一化
筆者還特地設計了一個 value t 可以讓我們選擇

畫出來的圖形跟原本的差不多
只是由於特徵縮放後 data_x 的範圍變了
線段的繪製範圍要稍作修改圖表才會比較好看

# 心得

最基礎的線性規劃大致上就是這樣
課程中其實還有教你如何顯示出 β₀和 β₁的修正過程
還有繪製動畫
甚至是推廣至高次多項式的回歸
但礙於版面因素筆者就不多作補充
身為小高二的筆者要應付一卡車的高三數學著實花了不少腦筋了
好在有先偷學一點微積分跟線性規劃 (汗
不然可能真的會應付不來呢

# Full Code

以下是筆者自己實作時的 code
由於功能有些許差異
可能與上面敘述有些出入就是了
(資料係課堂統一給定數據)

import numpy as np
import matplotlib.pyplot as plt
from IPython import display

data_x = np.array([24 ,22 ,15 ,4  ,9  ,20 ,5  ,3  ,17 ,19 ,13 ,10 ,12 ,11 ,16 ,27 ,16 ,16 ,6  ,20 ])
data_y = np.array([591,543,410,310,319,520,338,330,501,508,399,331,390,390,431,660,409,430,323,524])

def fc(x,t):
  return (x-np.min(x))/(np.max(x)-np.min(x)) if t==1 else (x-np.mean(x))/np.std(x)

def f(x):
  return b0 + b1 * x

def Loss(x, y):
    return 0.5 * np.sum((y - f(x)) ** 2)

t_x=fc(data_x[:len(data_x)//2],1)
t_y=data_y[:len(data_y)//2]
v_x=fc(data_x[len(data_x)//2:],1)
v_y=data_y[len(data_y)//2:]

b0=np.random.rand()
b1=np.random.rand()

ETA = 0.001
t_ct = 20000
LL=[]
Step0=[]
Step1=[]

def train():
  global ct
  global b0
  global b1
  for i in range(t_ct):
    b0 = b0 - ETA * np.sum((f(t_x) - t_y))
    b1 = b1 - ETA * np.sum((f(t_x) - t_y) * t_x)
    current_loss = Loss(t_x, t_y)
    ct += 1
    if not ct%500:
        LL.append(current_loss)
        Step0.append(b0)
        Step1.append(b1)

def check():
    x = np.linspace(0,1,2000)
    plt.figure(figsize=(8,6))
    plt.title("train-data",fontsize=20,color='white')
    plt.xticks(range(0,2),color='white')
    plt.yticks(range(0,2),color='white')
    plt.xlabel("x",fontsize=20,color='white')
    plt.ylabel("y",fontsize=20,color='white')
    plt.plot(x, f(x),label="mod func")
    plt.plot(t_x, t_y, 'go',label="train data")
    plt.plot(v_x, v_y, 'ro',label="valid data")
    plt.legend()
    plt.show()

def loss():
  plt.figure(figsize=(8,5))
  plt.title("change",fontsize=20,color='white')
  plt.xlabel("train times(x100)",fontsize=20,color='white')
  plt.ylabel("Loss",fontsize=20,color='white')
  plt.plot(LL,":o")
  plt.show()

def path():
  for i in range(len(Step0)):
    plt.figure(figsize=(8,5))
    plt.title("Movement path",fontsize=20,color='white')
    plt.xticks(range(0,2),color='white')
    plt.yticks(range(0,2),color='white')
    plt.xlabel("b0",fontsize=20,color='white')
    plt.ylabel("b1",fontsize=20,color='white')
    plt.scatter(Step0[:i], Step1[:i], s=20,color='r')
    plt.pause(0.1)
    display.clear_output(wait=True)

if __name__=="__main__":
  train()
  loss()
  path()
  check()