3.一進制線性回歸

線性回歸分析（Linear Regression Analysis)是确定兩種或兩種以上變量間互相依賴的定量關系的一種統計分析方法。對于一進制線性回歸而言，其模型主要假設為:

\[y=ax+b \tag{1}

\]

現在已經有了一組資料集\(\{(x^{(i)},y^{(i)})\}_{i=1}^{N}\),想要找到一個合适的參數\(a\)和\(b\) ，使得模型對每個\(x^{(i)}\)預測出來的值\(\hat{y}^{(i)}\)與真實值\(y^{(i)}\)盡可能的非常接近。基于這個目标，可以将均方誤差作為損失函數：

\[\begin{align}

L(w,b)

&=\frac{1}{2N}\sum_{i=1}^N(\hat{y}^{(i)}-y^{(i)})^2\\

&=\frac{1}{2N}\sum_{i=1}^N(ax^{(i)}+b-y^{(i)})^2

\end{align}

\tag{2}

\]

對于該損失函數，可以令其最小化求出合适的參數\(a\)和\(b\) 。這種基于均方誤差最小化求解線性回歸參數的方法稱為最小二乘法。從資料集到求出參數\(a\)和\(b\) 的過程稱為學習。所用到的資料集稱為知識積累。

\[\begin{align}

a^*,b^* &= \arg \max_{a,b}\mathrm{L}(w,b)

\end{align}

\tag{3}

\]

對于以上優化問題，可以求相應的偏導數令其為\(0\)得到相應的解析解，也可以通過梯度下降法去求解，下面為求解過程：

\[\begin{align}

\frac{\partial\mathrm{L}(w,b)}{\partial a}

&=\frac{1}{N}\sum_{i=1}^N\{(ax^{(i)}+b-y^{(i)})x^{(i)}\} \\

&=\frac{1}{N}\sum_{i=1}^N\{({ax^{(i)}}^2+bx^{(i)}-y^{(i)}x^{(i)})\} \\

&=\frac{1}{N} \{ a\sum_{i=1}^N{x^{(i)}}^2+b\sum_{i=1}^Nx^{(i)}- \sum_{i=1}^Ny^{(i)}x^{(i)} \}

\end{align}

\tag{4}

\]

\[\begin{align}

\frac{\partial\mathrm{L}(w,b)}{\partial b}

&=\frac{1}{N}\sum_{i=1}^N(ax^{(i)}+b-y^{(i)}) \\

&=\frac{1}{N} \{ a\sum_{i=1}^N{x^{(i)}}+Nb- \sum_{i=1}^Ny^{(i)} \}

\end{align}

\tag{5}

\]

對于(4)和(5),分别令其為\(0\),可得到\(w^*\)和\(b^*\)：由于求解過程不同，可能得到的表達式也會不同

\[w^*=\frac{\sum_{i=1}^N(x^{(i)}-\bar{x})(y^{(i)}-\bar{y})}{\sum_{i=1}^N(x^{(i)}-\bar{x})^2}

\tag{6}

\]

\[b^{*}=\bar{y}-a\bar{x} \tag{7}

\]

通過觀察可以發現，\(w^*\)的分子是自變量和因變量的協方差，分母為自變量的方差。如果對資料做中心化處理，那麼可以簡化以上的結果：

\[w^*=\frac{\sum_{i=1}^Nx^{(i)}y^{(i)}}{\sum_{i=1}^N{x^{(i)}}^2}

\tag{8}

\]

\[b^*=0 \tag{9}

\]

下面通過代碼示範基于梯度下降法來求解的過程：

def get_data(true_a,true,b):
    x=np.linspace(0,2,100)
    y=a*x+b+np.random.normal(size=len(x))
    return(x,y)
def lin_model(x,a,b):
    return a*x+b
def loss(X,Y,a,b):
    return np.sum(Y-(a*X+b))/(2*len(X))
def loss_da(X,Y,a,b):
    return ((a*np.sum(X**2)+b*np.sum(X)-np.sum(X*Y)))/len(X)
def loss_db(X,Y,a,b):
    return (a*np.mean(X)+b-np.mean(Y))/len(X)
def SGD(X,Y,a_init,b_init,lr=0.01,error=1e-10,maxIter=100000):
    a,b=(a_init,b_init)
    for iter in range(0,maxIter):
        da=loss_da(X,Y,a,b)
        db=loss_db(X,Y,a,b)
        new_a=a-lr*da
        new_b=b-lr*db
        if(np.abs(loss(X,Y,a,b)-loss(X,Y,a_init,b_init))<error):
            return (new_a,new_b)
        else:
            a,b=(new_a,new_b)
    return (new_a,new_b)
           

true_a,true_b=(2,4)
X,Y=get_data(true_a,true_b)
SGDLinerModel(X,Y,1,2)
#(1.9253542219938329, 4.022734902793057)