不得不說安卓老師是真的用心良苦，給我們把程式設計題弄成了填空題，但是很可惜原版使用的是Octave和MATLAB所以作為初學者我就直接當程式設計題用python去做了

問題：讓你繪製一個5階單位陣

答案：

在本練習的這一部分，你將應用線性迴歸來實現通過用一個變數x（人口）預測y(食品卡車的利潤)。假設你是一位執行長，餐飲連鎖店正在考慮不同的城市開設新的出口。這個連鎖已經把卡車開到了各個城市，現在你有來自城市的利潤和人口的資料了。你希望使用這些資料來幫助你選擇擴充套件哪個城市。

你希望通過這些資料來輔助你決定下一步選擇擴充套件哪個城市。

檔案EX1DATA1.TXT

在開始任何任務之前，通過資料視覺化來理解資料是有很有用的。對於這個資料集，可以使用散點圖來視覺化他。因為它只只有兩個屬性需要繪製（利潤和人口）。

而在現實生活中你會遇到很多其他問題是不能在二維圖上繪製的。

在EX1.M中，將資料集從資料檔案載入到變數X中和Y：

問題：讓你繪製一個數據集的散點圖

答案：

在這一部分中，你將採用梯度下降法通過線性迴歸找到合適的引數θ去擬合到我們的資料集。

這裡題目寫的是θ^tx就是 1*2  *  2*m 我們為了程式設計中方便觀察（如果

為資料新增θ0這個截距項（線性方程的常數項）的維度並將θ初始化為0，設定學習速率為0.01，迭代次數為1500

當你執行梯度下降最小化代價函式時，通過計算代價值將幫你監控收斂性，計算代價函式 在你做的時候切記x和y不是標量而是代表訓練集的行矩陣，當θ為初始值時這個函式執行出來的結果是32.07

接下來你將對文字執行梯度下降 安卓老師已經幫我們寫好了迴圈結構 只需要我們填空（再次為老師點贊是真的良心）我們想要把J(θ)降到最小 當你更新引數時切記是更新x、y的向量而不是一個x、y的引數值 如果你忘了怎麼做

吳恩達再次強調你要做的是向量運算而不是數值運算

為了更好地理解代價函式 現在你將繪製θ0和θ1二維網格圖。在這部分 你不需要編寫任何新的程式碼 但是你應該理解你是如何這些程式碼是如何繪製這些影象的。

在程式碼執行後會生成上述影象

這些圖表的目的是向你展示J (θ)是如何隨θ0和θ1的變化而變化的。成本函式J (θ)是碗形的 具有全域性最小值(在等高線圖中比在3D曲面圖更容易看出)這個最小值對應的是θ0和θ1就是θ0、θ1的最優點 並且梯度下降的每一步都更接近這個點。

完整程式碼及註解

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def warmupExercise():
    E5 = np.eye(5)  # eye(5)代表5階單位陣
    print('這是一個五階單位陣')
    print(E5)

def plotDatainit(x, y):

    plt.title('ex1Sample')#圖的標題
    plt.plot(x, y, 'rx', ms=5)#ms=5代表散點大小為5個單位
    plt.xticks(np.arange(4,24,2))#設定x軸的尺度（這裡和吳老師在pdf上繪製一致從4開始24結束，2為組距）
    plt.yticks(np.arange(-5,25,5))#設定y軸尺度
    plt.xlabel('Population of City in 1,0000s')#橫軸的含義：x:城市人口
    plt.ylabel('Profit in $1,0000s')#縱軸含義y:收益
    plt.show()

def plotData(x, y,theta0,theta1):

    plt.title('ex1Sample')#圖的標題
    plt.plot(x, y, 'rx', ms=5)#ms=5代表散點大小為5個單位
    plt.xticks(np.arange(4,26,2))#設定x軸的尺度（這裡和吳老師在pdf上繪製一致從4開始24結束，2為組距）
    plt.yticks(np.arange(-5,25,5))#設定y軸尺度
    plt.xlabel('Population of City in 1,0000s')#橫軸的含義：x:城市人口
    plt.ylabel('Profit in $1,0000s')#縱軸含義y:收益
    x = np.arange(4, 24)
    y = theta0+theta1*x
    plt.plot(x, y)
    plt.show()

def computeCost(X,y,theta,num):
  m=num
  result=np.sum((np.dot(X,theta)-y)**2)/(2*m)#計算假設函式的值
  return result

def gradientDescent(X,y,theta,num,studyrate,iterations):
    m=num
    J_every=np.zeros((iterations,1))#用於儲存每次迭代計算出的costfunction值
    theta_every=np.zeros((interations,2))
    XT=X.T#為保證矩陣乘法行列要求 接下來計算會用到X的轉置
    for i in range(interations):
        #dJdivdtheta = np.sum(np.dot(XT,((np.dot(X, theta) - y))/m )) #這是錯的
        dJdivdtheta = (np.dot(XT,((np.dot(X, theta) - y))/m ))#這是對的
        theta = theta - studyrate * dJdivdtheta
        theta_every[i][0]=theta[0][0]
        theta_every[i][1]=theta[1][0]
        J_every[i] = computeCost(X, y, theta,num)
    return theta,J_every,theta_every

def showCostfunction(interations):
  for i in range(interations):
    print('第',i+1,'次迭代的代價函式值為',costhistory[i],'theta0和theta1分別是',thetahistory[i])

def predictProfit():
  predict1=np.dot(np.array([1,3.5]),thetai)#預測3.5W人口的利潤
  predict2=np.dot(np.array([1,7]),thetai)#預測7W人口的利潤
  print("3W5人口利潤是",(predict1[0]*10000))
  print("7W人口的利潤是",(predict2[0]*10000))

def plot3DJtheta():
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    xx=list(thetahistory[:,0])#把thetahistory的第1列即使theta0變成一個1500個元素組成的列表
    yy=list(thetahistory[:,1])#把thetahistory的第2列即使theta1變成一個1500個元素組成的列表
    zz=costhistory.reshape(1500)#把costhistory形狀弄成1500個的橫著的
    # 網上抄的如上
    ax.set_xlabel("theta0", color='b')
    ax.set_ylabel("theta1", color='g')
    ax.set_zlabel("Z(theta)", color='r')
    ax.plot_trisurf(xx, yy, zz,color='red')
    plt.show()

def plotContour(thetai,X,y,num):#這個函式是網上抄的-0-
    theta0_vals = np.linspace(-10, 10, 100)
    theta1_vals = np.linspace(-1, 4, 100)
    J_vals = np.zeros((len(theta0_vals), len(theta1_vals)))
    for i in range(len(theta0_vals)):
        for j in range(len(theta1_vals)):
            t = np.array([theta0_vals[i], theta1_vals[j]]).reshape((2, 1))
            J_vals[i, j] = computeCost(X, y, t, num)
    J_vals = J_vals.T
    # fig = plt.figure()
    # ax = Axes3D(fig)
    # ax.plot_surface(theta0_vals, theta1_vals, J_vals, rstride=1, cstride=1, cmap='rainbow')
    # plt.xlabel('theta_0')
    # plt.ylabel('theta_1')
    # 這部分是3D圖
    # 填充顏色，20是等高線分為幾部分
    plt.contourf(theta0_vals, theta1_vals, J_vals, 20, alpha=0.6, cmap=plt.cm.hot)
    plt.contour(theta0_vals, theta1_vals, J_vals, colors='black')
    plt.plot(thetai[0], thetai[1], 'r', marker='x', markerSize=10, LineWidth=2)  # 畫點
    plt.show()


warmupExercise()
dataset=np.loadtxt('ex1data1.txt',delimiter=',')#讀取資料集文字
xx=dataset[:,0]#取出文字的第一列
yy=dataset[:,1]#取出文字的第二列
num=len(xx)
#print(num) #看看有多少個x (有97個數據)
x=xx.reshape(num,1)#轉成num行的列向量
y=yy.reshape(num,1)#同上
thetainit=np.zeros((2,1))#由題設theta包含兩個引數常數項theta0和一階係數theta1 後面會更新theta所以初始化是一個空矩陣 是2行1列的矩陣
juzhenx=np.array(x).reshape(len(x),1)#資料集中屬性x轉換成一個n行1列的矩陣
X0=np.ones((len(x), 1))#因為假設函式是theta0加theta1*x 所以我們後續的計算涉及到矩陣運算
X = np.c_[X0, juzhenx]#把矩陣x(97行1列)拼在X0(97行1列)右邊 形成一個97行2列的矩陣X
#X=np.hstack((X0,juzhenx))#這是矩陣拼接的另外一個方法
#print(X)#看看X是否符合要求
#print(X.shape)#看看X行列是否符合要求
studyrate=0.01
interations=1500
print("Costfuction的初始值是",computeCost(X,y,thetainit,num))#檢視還未更新的代價函式值
thetai,costhistory,thetahistory=(gradientDescent(X,y,thetainit,num,studyrate,interations))
showCostfunction(interations)
predictProfit()
plotData(x, y,thetai[0],thetai[1])
plot3DJtheta()
plotContour(thetai,X,y,num)

見下