前言：學(xué)習(xí)了吳恩達(dá)機(jī)器學(xué)習(xí)線性回歸部分內(nèi)容后，再結(jié)合很多優(yōu)秀博客總結(jié)，搬運(yùn)出來的干貨。

1、原理

2、單變量線性回歸

數(shù)據(jù)準(zhǔn)備：https://blog.csdn.net/Carl_changxin/article/details/100824723，第一列當(dāng)作人口Population，第二列當(dāng)作收益Profite

實(shí)現(xiàn)技術(shù)：梯度下降

實(shí)現(xiàn)目標(biāo)：

（1）作原始數(shù)據(jù)散點(diǎn)圖

（2）作線性回歸模型折線圖

（3）根據(jù)模型預(yù)測(cè)一些數(shù)據(jù)

（4）作代價(jià)函數(shù)隨迭代次數(shù)增加而變化的折線圖；會(huì)發(fā)現(xiàn)代價(jià)函數(shù)慢慢趨向某個(gè)最小值，從而調(diào)整迭代次數(shù)

（5）作代價(jià)函數(shù)和參數(shù)的3d圖；會(huì)發(fā)現(xiàn)代價(jià)函數(shù)屬于凸函數(shù)，只有全局最小值

（6）作代價(jià)函數(shù)和參數(shù)的等高線圖；即由3d圖投影所得，可觀察參數(shù)變化過程

實(shí)現(xiàn)代碼：如下

            
              # 1、單變量線性回歸python實(shí)現(xiàn)
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import LogNorm
from mpl_toolkits.mplot3d import axes3d, Axes3D

'''
分析難點(diǎn)：
（1）hypothesis = X * theta.T
X = [x0 x1 x2 ... xn]，其中x0~xn是m*1的列向量，所有X是m*(n+1)的矩陣；
X是包含所有訓(xùn)練數(shù)據(jù)的特征矩陣，每一行都代表一條數(shù)據(jù)
theta = [theta0 theta1 theta2 ... theta-n]，1*(n+1)參數(shù)向量
所以hypothesis的值是m*1的向量，包含所有測(cè)試值的向量
（2）theta = theta - (alpha/y.size) * ((X * theta.T - y).T * X)
根據(jù)上面的分析X * theta.T是m*1的向量，包含所有測(cè)試值的向量；
(X * theta.T - y).T就是1*m的行向量，X是m*(n+1)的矩陣，故相乘結(jié)果為1*(n+1)的向量；
'''


# 代價(jià)函數(shù);X=[x0 x1 ...xn]是m*(n+1)的矩陣,y=[m*1列向量]，theta=[1*(n+1)行向量]
def computeCost(X, y, theta):
    # 假設(shè)函數(shù);由于X中包含所有訓(xùn)練數(shù)據(jù)，所以得到的結(jié)果是所有測(cè)試結(jié)果組成的m*1列向量
    hypothesis = X * theta.T
    # 訓(xùn)練數(shù)據(jù)個(gè)數(shù),或者用m = y.shape[1]
    m = y.size
    # 代價(jià)函數(shù)
    J = 1 / (2 * m) * np.sum(np.power(hypothesis - y, 2))
    return J


# 梯度下降函數(shù);X=[x0 x1 ...xn]是m*(n+1)的矩陣,y=[m*1列向量]，theta=[1*(n+1)行向量]
def gradientDescent(X, y, theta, alpha, iters):
    # cost用來記錄迭代每一次的代價(jià)函數(shù)值；用長度為iters的數(shù)組記錄；初始化為0
    cost = np.zeros(iters)
    # thetaRecord記錄每一次迭代的theta;thetaRecord是二維數(shù)組
    thetaRecord = np.zeros([iters, 2])
    # 每迭代一次，就要循環(huán)更新一次所有參數(shù)的值
    for i in range(iters):
        thetaRecord[i] = np.array(theta)[0]
        theta = theta - (alpha / y.size) * ((X * theta.T - y).T * X)  # 更新theta向量
        cost[i] = computeCost(X, y, theta)
    return theta, cost, thetaRecord


if __name__ == '__main__':
    # ---繪圖設(shè)置參數(shù)---
    font = {
        'family': 'SimHei',
        'weight': 'bold',
        'size': '16'
    }
    plt.rc('font', **font)
    plt.rc('axes', unicode_minus=False)

    # （1）讀取訓(xùn)練數(shù)據(jù)集，并繪制數(shù)據(jù)集散點(diǎn)圖
    train_data = pd.read_csv('ex1data1.txt', names=['Population', 'Profit'])
    train_data.plot(kind='scatter', x='Population', y='Profit', figsize=(12, 8))
    plt.show()

    # 初始化變量
    train_data.insert(0, 'Ones', 1)
    X = train_data.iloc[:, [0, 1]]
    y = train_data.iloc[:, 2]
    X = np.mat(X.values)  # X是矩陣
    y = np.mat(y.values)
    y = y.T  # 將y變成列向量
    theta = np.mat('0 0')
    iters = 5000
    alpha = 0.01
    bestPara, cost, thetaRecord = gradientDescent(X, y, theta, alpha, iters)

    # （2）畫出線性回歸函數(shù)
    x = np.linspace(train_data.Population.min(), train_data.Population.max(), 100)
    f = bestPara[0, 0] + (bestPara[0, 1] * x)
    fig, ax = plt.subplots(figsize=(12, 8))
    # 將兩個(gè)圖像繪制在同一個(gè)figure上
    ax.plot(x, f, 'r', label='預(yù)測(cè)函數(shù)')  # 預(yù)測(cè)函數(shù)折線圖
    ax.scatter(train_data.Population, train_data.Profit, label='訓(xùn)練數(shù)據(jù)')  # 訓(xùn)練數(shù)據(jù)散點(diǎn)圖
    ax.legend(loc=2)
    ax.set_xlabel('人口')
    ax.set_ylabel('利潤')
    ax.set_title('預(yù)測(cè)利潤和人口數(shù)量')
    plt.show()

    # （3）預(yù)測(cè)人口規(guī)模為3.5萬和7萬的利潤值
    predict1 = bestPara[0, 0] * 1 + (bestPara[0, 1] * 3.5)
    print('當(dāng)人口為35,000時(shí),我們預(yù)測(cè)利潤為', predict1 * 10000)
    predict2 = bestPara[0, 0] * 1 + (bestPara[0, 1] * 7)
    print('當(dāng)人口為70,000時(shí),我們預(yù)測(cè)利潤為', predict2 * 10000)

    # （4）畫出代價(jià)函數(shù)隨迭代次數(shù)增長而變化的圖像
    # subplots的返回值有figure, Axes對(duì)象
    fig, ax = plt.subplots(figsize=(12, 8))
    ax.plot(np.arange(iters), cost, 'r')  # 繪圖
    ax.set_xlabel('迭代次數(shù)')
    ax.set_ylabel('輸出代價(jià)')
    ax.set_title('誤差和訓(xùn)練狀態(tài)')
    plt.show()

    # （5）畫出代價(jià)函數(shù)和參數(shù)的3d圖像和等高線圖
    theta0_vals = np.linspace(-10, 10, 100)
    theta1_vals = np.linspace(-1, 4, 100)
    xs, ys = np.meshgrid(theta0_vals, theta1_vals)
    J_vals = np.zeros(xs.shape)
    for i in range(0, theta0_vals.size):
        for j in range(0, theta1_vals.size):
            t = np.array([theta0_vals[i], theta1_vals[j]])
            t = np.mat(t)
            J_vals[i][j] = computeCost(X, y, t)
    J_vals = np.transpose(J_vals)  # 轉(zhuǎn)置函數(shù)，轉(zhuǎn)成行向量
    # 繪制3D圖
    fig1 = plt.figure(2)
    ax = fig1.gca(projection='3d')  # 將figure變?yōu)?d
    ax.plot_surface(xs, ys, J_vals)  # 繪制3D圖
    plt.xlabel(r'$\theta_0$')
    plt.ylabel(r'$\theta_1$')
    # 繪制等高線,觀察迭代過程中參數(shù)的變化
    plt.figure(3)
    lvls = np.logspace(-2, 3, 20)  # 表示20**-2,20**-1,...20**3
    plt.contour(xs, ys, J_vals, levels=lvls, norm=LogNorm())  # 繪制等高線
    # for i in range(iters):
    #     plt.plot(thetaRecord[i][0], thetaRecord[i][1], c='r', marker="x")

    # 若不想觀察所有的點(diǎn)，可以每隔500次迭代繪制一次
    for i in range(0, iters, 500):
        plt.plot(thetaRecord[i][0], thetaRecord[i][1], c='r', marker="x")
    # 最后一次，手動(dòng)輸出
    plt.plot(thetaRecord[iters - 1][0], thetaRecord[iters - 1][1], c='r', marker="x")
    plt.show()

            
          

實(shí)現(xiàn)截圖：

?3、多元線性回歸

數(shù)據(jù)準(zhǔn)備：https://blog.csdn.net/Carl_changxin/article/details/100824723，

實(shí)現(xiàn)技術(shù)：梯度下降，特征縮放（均值歸一化）

實(shí)現(xiàn)目標(biāo)：

（1）訓(xùn)練數(shù)據(jù)，得到最佳參數(shù)theta；測(cè)試數(shù)據(jù)，預(yù)測(cè)結(jié)果

實(shí)現(xiàn)代碼：

            
              # 1、二元線性回歸python實(shí)現(xiàn)
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import LogNorm
from mpl_toolkits.mplot3d import axes3d, Axes3D


# 代價(jià)函數(shù);X=[x0 x1 ...xn]是m*(n+1)的矩陣,y=[m*1列向量]，theta=[1*(n+1)行向量]
def computeCost(X, y, theta):
    # 假設(shè)函數(shù);由于X中包含所有訓(xùn)練數(shù)據(jù)，所以得到的結(jié)果是所有測(cè)試結(jié)果組成的m*1列向量
    hypothesis = X * theta.T
    # 訓(xùn)練數(shù)據(jù)個(gè)數(shù),或者用m = y.shape[1]
    m = y.size
    # 代價(jià)函數(shù)
    J = 1 / (2 * m) * np.sum(np.power(hypothesis - y, 2))
    return J


# 梯度下降函數(shù);X=[x0 x1 ...xn]是m*(n+1)的矩陣,y=[m*1列向量]，theta=[1*(n+1)行向量]
def gradientDescent(X, y, theta, alpha, iters):
    # cost用來記錄迭代每一次的代價(jià)函數(shù)值；用長度為iters的數(shù)組記錄；初始化為0
    # cost = np.zeros(iters)
    # 每迭代一次，就要循環(huán)更新一次所有參數(shù)的值
    for i in range(iters):
        theta = theta - (alpha / y.size) * ((X * theta.T - y).T * X)  # 更新theta向量
        # cost[i] = computeCost(X, y, theta)
    return theta


# 特征縮放
def featureNormalize(X):
    X_norm = X
    mu = np.zeros((1, X.shape[1]))
    sigma = np.zeros((1, X.shape[1]))
    for i in range(X.shape[1]):
        mu[0, i] = np.mean(X[:, i])  # 均值
        sigma[0, i] = np.std(X[:, i])  # 標(biāo)準(zhǔn)差
    X_norm = (X - mu) / sigma
    return X_norm, mu, sigma


if __name__ == '__main__':
    # 讀取訓(xùn)練數(shù)據(jù)集
    train_data = pd.read_csv('ex1data2.txt', names=['id', 'feature1', 'feature2', 'feature3', 'label'])
    # 刪除id列
    train_data = train_data.drop(columns=['id'])

    # 初始化變量
    X = train_data.iloc[:, [0, 1, 2]]
    y = train_data.iloc[:, 3]
    X = np.mat(X.values)  # X從DataFrame類型轉(zhuǎn)化成矩陣
    y = np.mat(y.values)
    y = y.T  # 將y變成列向量
    # 特征縮放
    X, mu, sigma = featureNormalize(X)
    # 特征縮放后再新增x0列
    b = np.ones(y.size)
    X = np.insert(X, 0, values=b, axis=1)
    theta = np.mat(np.zeros(X.shape[1]))
    iters = 5000
    alpha = 0.01

    # (1)訓(xùn)練數(shù)據(jù)，得到最佳參數(shù)theta；測(cè)試數(shù)據(jù)，預(yù)測(cè)結(jié)果
    # 訓(xùn)練
    bestPara = gradientDescent(X, y, theta, alpha, iters)
    # 測(cè)試
    data = [151.5, 41.3, 58.5]
    testx = np.mat(data)
    # 對(duì)測(cè)試數(shù)據(jù)進(jìn)行特征縮放
    testx = ((testx - mu) / sigma)
    # 對(duì)測(cè)試數(shù)據(jù)添加x0列
    b = np.ones(testx.shape[0])
    testx = np.insert(testx, 0, values=b, axis=1)
    # print(testx)
    # print(bestPara)
    # dot方法，返回兩個(gè)數(shù)組的內(nèi)積;返回矩陣的乘積
    price = testx.dot(bestPara.T)
    # predit value is 17.597830
    print('predit value is %f ' % (price))



            
          

?

?

?

?

?

?

?