【參考資料】 【1】《概率論與數理統計》 【2】 http://scikit-learn.org /stable/auto_examples/ linear_model/ plot_ols.html # sphx-glr-auto-examples-linear-model-plot-ols-py

數學基礎

1. 基本定義： 隨機變數Y（因變數）和變數X（自變數）之間存在某種關係。假設隨機變數Y的期望存在且是X的函式$\mu=\mu(x)$,這裡$\mu(x)$稱為Y關於X的迴歸函式，簡稱Y關於X的迴歸。

如上圖所表示，如果$\mu=\mu(x)$

$\begin{cases} y = a + bx + \varepsilon \\ \varepsilon \sim N(0, \sigma^2) \end{cases}$ 備註：相當於一個線性的模型上疊加了一個方差很小的正態分佈誤差

2. 推導過程

已知$Y_i\sim N(a+bX_i, \sigma^2)$,由於$Y_i$獨立，因此根據極大似然估計法，構建聯合分佈密度函式如下：

$L$

上式取最大值，等價於$Q(a,b)=(y_i - a - bx_i)^2$取最小值，因此對其求偏導數

$\begin{cases} \dfrac{\partial Q}{\partial a} = -2 \sum\limits_{i=1}^n(y_i-a-bx_i)=0 \\ \dfrac{\partial Q}{\partial b} = -2 \sum\limits_{i=1}^n(y_i-a-bx_i)x_i=0 \end{cases}$

程式實現 (基於sklearn)

# -*- coding: utf-8 -*-
import numpy  as np
import matplotlib.pyplot as plt
from sklearn import linear_model


def _test_one_lr():

    index  = np.linspace(0, 10, 20);

    x_test = np.hstack((index, np.ones(20))).reshape(2,20)

    """
    隨機變數為y = 0.3x + 0.2 + delta
    從矩陣的角度看y = wx
    w = [0.3,0.2]^T x = [x, 1]^T
    """
    delta  = np.random.random_sample(20)

    y_test = 0.3*index + 0.2*np.ones(20) + delta

    regr = linear_model.LinearRegression()

    """
    輸入時轉置，由行矩陣變成列矩陣
    fit方法為進行模型訓練
    """
    regr.fit(np.transpose(x_test),np.transpose(y_test))

    """
    根據訓練得到的模型進行預測
    """
    y_pred = regr.predict(np.transpose(x_test))

    plt.scatter(index, y_test,  color='black')
    plt.plot(index, y_pred, color='blue', linewidth=2)

    plt.xticks(())
    plt.yticks(())

    plt.show()


if __name__ == "__main__":

    _test_one_lr()

輸出結果：

程式實現 (使用極大似然估計公式推導)

# -*- coding: utf-8 -*-
import numpy  as np
import matplotlib.pyplot as plt 



def _test_one_lr(): 

    index  = np.linspace(0, 10, 20);

    """
    隨機變數為y = 0.3x + 0.2 + delta
    從矩陣的角度看y = wx
    w = [0.3, 0.2]^T x = [x, 1]^T

    delta 是一個混淆變數
    """
    delta  = np.random.random_sample(20)/10

    y_test = 0.3*index + 0.2*np.ones(20) + delta


    x_mean = np.mean(index)
    y_mean = np.mean(y_test)

    b_up   = np.dot(index, y_test) - 20 * x_mean * y_mean

    b_down = np.dot(index, index) - 20 * x_mean * x_mean

    b = b_up / b_down

    a = y_mean - (b * x_mean)

    
    print("function: y =  (%f) * x + %f" %(b, a))

    #輸出: function: y =  (0.300204) * x + 0.252471

"""
說明：

一元線性迴歸的最小二乘法實現，採用極大似然估計的概率公式。對應的筆記《最小二乘法（一元線性迴歸）》

作者：fredric

日期：2018-8-11

"""
if __name__ == "__main__":

    _test_one_lr()