基于矩陣向量的單變數線性回歸（python實作 ）

梯度下降法實作

在本部分的練習中，我們將使用一個變數實作線性回歸，我們需要找到城市人口數量與城市餐飲利潤之間的關系，然后根據城市人口數量來預測城市餐飲的利潤，資料集
鏈接：https://pan.baidu.com/s/1Zkdtsgvan8-YhwATFcBbSw
提取碼：ABi8

匯入用到的模塊以及資料，進行可視化

#匯入用到的模塊
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d #mpl_toolkits.mplot3d是matplotlib中提供的專門畫三維圖的工具包
from mpl_toolkits.mplot3d import Axes3D #匯入Axes3D用以創建三維圖

#匯入資料，進行資料的可視化
data = pd.read_csv('ex1data1.txt',names=['population','profit'])#當資料與py檔案在同一目錄下，可以直接用檔案名匯入資料
data.plot.scatter( 'population','profit',c='r',label='population')#可以使用列名進行資料讀取
plt.title('Training Set')
plt.show()

資料可視化圖

代價函式

(1)假設函式：H_θ（x）=θ₀+θ₁ ? x₁
我們增加一個虛擬的輸入特征x₀ = （x₀₁ = 1，x₀₂ = 1，…x_0n = 1）， n∈N,此時H_θ（x）= θ^Tx
我們的目的是求出θ0，θ1的值，使得代價函式可以取最小值，
(2) 代價函式

代價函式

將其轉化為矩陣向量的形式 =

#定義代價函式
def costFunction(X,y,theta):
    inner = np.power(X*theta.T - y,2)#np.power(a,3)就是a的立方
    return np.sum(inner)/(2*len(X))

梯度下降

梯度下降是一種用來求代價函式最小值的一種演算法，

重復計算這個公式，直到函式收斂，（注意更新值θ應同時更新）
其中的α稱之為步長，在最優化課中我們可以調整α的值，
這個α如果過小，則收斂很慢； 如果過大，則可能導致不收斂，

對于J(θ)的偏導數學推導程序如下：

其中x是向量，x_j是特征向量x中的第j元素，

#定義梯度下降函式
def gradientDescent(X,y,theta,alpha,iters):
    costs =[]
    for i in range(iters):
        wucha= (X*theta.T - y)
        theta -= (wucha.T*X)*alpha/len(X)  # theta.T代表矩陣的轉置
        cost = costFunction(X,y,theta)
        costs.append(cost)
##      if i %100==0:   沒迭代一百次輸出一個cost
##         print(cost)
    return  theta,costs

完整代碼實作

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d #mpl_toolkits.mplot3d是matplotlib中提供的專門畫三維圖的工具包
from mpl_toolkits.mplot3d import Axes3D #匯入Axes3D用以創建三維圖

# 匯入資料pandas
data = pd.read_csv('ex1data1.txt',names=['population','profit'])#names是添加列名
data.plot.scatter( 'population','profit',c='r',label='population')
plt.title('Training Set')
plt.show()

#資料處理
data.insert(0,'ones',1)#對x要插入一列，列名為ones，值全部是一
X = data.iloc[:,0:2]#0列和1列所有元素，不包括第二列，或者2變為-1，代表最后一列（有兩個冒號），也取不到
##cols=data.shape[1]  顯示data的列數，.[0]是行數,或者用這個cols-1來切片
y = data.iloc[:,2:]#或者[:,-1]只有一個冒號就是所有行最后一列，因為第二列就是最后一列，所以：后邊可以省略一個數字但是此時需要.reshape(9)
##print(len(X))輸出X的樣本數，使用[:,-1]是一維元組，此時需要y=y.values然后y=y.reshape(97,1)改為二維

##當使用-1的時候
##y = data.iloc[:,-1]
##y=y.values
##y=y.reshape(97,1)
##print(y.shape)


#將資料X,Y的資料結構（dataframe）轉換為陣列（ndarray）
X = np.matrix(X.values) #把X與Y轉換為矩陣，方便接下來的矩陣操作X.values轉換為array的形式，
y = np.matrix(y.values) #matrix的優勢就是相對簡單的運算子號，比如兩個矩陣相乘，就是用符號*，但是array相乘不能這么用，得用方法.dot()
#print(y.shape) 檢查X,y的維度


#定義代價函式
def costFunction(X,y,theta):
    inner = np.power(X*theta.T - y,2)#np.power(a,3)就是a的立方
    return np.sum(inner)/(2*len(X))

theta =np.zeros((1,2))#或者用這個theta = np.matrix(np.array([0,0]))，都是生成1*2矩陣(代價函式theta有轉置)，若沒有轉置使用np.zeros((2,1))
##cost_init=costFunction(X,y,theta)
##print(cost_init)    輸出初始的代價函式值


#定義梯度下降函式
def gradientDescent(X,y,theta,alpha,iters):
    costs =[]
   
    for i in range(iters):
        wucha= (X*theta.T - y)
        theta -= (wucha.T*X)*alpha/len(X)
        cost = costFunction(X,y,theta)
        costs.append(cost)
##      if i %100==0:   沒迭代一百次輸出一個cost
####        print(cost)
    return  theta,costs


alpha = 0.02
iters=2000
theta,costs=gradientDescent(X,y,theta,alpha,iters)

#  畫線性擬合后的圖形
##plt.figure() #創建畫布
##plt.subplot(1,2,1)   #一行輸出兩個圖形,且在第一個位置，畫另一個圖形時（1,2,2）則是同時輸出兩張圖
x=np.linspace(X[:,-1].min(), X[:,-1].max(), 100) #抽X向量中100個樣本，取樣本可以直接切片X[:,1:]=X[:,-1]都是X最后一列
f=theta[0,0] + (theta[0,1]*x) 
plt.plot(x,f,'r',label='Prediction',)
plt.scatter(X[:,1:2].tolist(),y.tolist(),label='Training data')  #畫散點圖的時候，X,和y都是多維的，用.tolist改成一維
plt.xlabel('Population')
plt.ylabel('Profit')
plt.legend(loc=2)#顯示標簽位置
plt.title('Predicted Profit vs. Population Size')
plt.show()

#輸出迭代次數與代價函式的變化
plt.plot(np.arange(iters),costs,label='cost vs iters')
plt.xlabel('iters')
plt.ylabel('cost')
plt.title('cost vs iters')
plt.legend(loc=1)
plt.show()
print(theta,costs[-1])


# 代價函式三維圖
def CostContour(X,y):
    J_vals = np.zeros((100,100))
    theta0_vals = np.linspace(-10,10,100)
    theta1_vals = np.linspace(-1,4,100)
    for i in range(len(theta0_vals)):
        for j in  range(len(theta1_vals)):
            z = np.matrix((theta0_vals[i],theta1_vals[j]))
            J_vals[i,j] = costFunction(X,y,z)
    return theta0_vals,theta1_vals,J_vals


#繪制三維圖        
fig=plt.figure()	
ax=Axes3D(fig)   #為創建的圖片視窗添加一個型別為3D的圖片，命名為ax
a,b,c = CostContour(X,y)
ax.plot_surface(a,b,c,cmap="rainbow")#cmap="rainbow"代表彩虹色
ax.set_xlabel('theta_0')
ax.set_ylabel('theta_1')
ax.set_zlabel('J(w,b)')
plt.show()

#繪制等高線圖
plt.figure() #畫等高線
lvls = np.logspace(-2, 3, 20) #以10**-2次方（0.01）和10**3（1000）之間的對數間隔繪制20條等高線，
a,_,c = CostContour(X,y)
_,b,_ = CostContour(X,y)
z = plt.contour(a,b,c,levels=lvls)
plt.clabel(z)  #顯示每條等高線的代價函式值
plt.show() #顯示影像

運行結果如下圖所示

離散資料可視化

線性擬合圖形

代價函式三維圖

等高線圖