# -*- coding: utf-8 -*-
"""
Created on Mon Jun 29 16:01:13 2020
@author: cheetah023
"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
#函数定义
#输出对角矩阵
def warmUpExercise(n):
A = np.eye(n)
return A
#损失函数
def computeCost(X, y, theta):
J = np.sum((X.dot(theta) - y) ** 2) / (2 * len(y))
return J
#梯度下降函数
def gradientDescent(X, y, theta, alpha, num_iters):
m = X.shape[0]
J_history = np.zeros((num_iters,1))
for i in range(1,num_iters):
h = np.dot(X,theta)
theta = theta - (alpha / m) * np.dot(X.T,h-y)
J_history[i] = computeCost(X, y, theta)
return theta,J_history
#Part 1: Basic Function
print('Part1.\n')
A = warmUpExercise(5)
print(A)
#Part 2: Plotting
print('Part2.\n')
data = np.loadtxt('ex1data1.txt',delimiter=',')
X = data[:,0]
y = data[:,1:2]#使用data[:,1]的话，y.shape打印出来是(97,)
print('y:')
print(y.shape)
plt.scatter(X, y,marker='x',c='r',)
plt.ylabel('Profit in $10,000s')
plt.xlabel('Population of City in 10,000s')
#Part 3: Cost and Gradient descent
print('Part3.\n')
m = len(y)
ones = np.ones(m)
theta = np.zeros((2,1))
X = np.column_stack((ones,X))#在X的首列插入全是1的列
print('X:')
print(X.shape)
iterations = 1500
alpha = 0.01
J = computeCost(X, y, theta)
print('With theta = [0 ; 0],Cost computed = %f' %(J))
print('Expected cost value (approx) 32.07\n')
J = computeCost(X, y, np.array([[-1],[2]]))
print('With theta = [-1 ; 2],Cost computed = %f' %(J))
print('Expected cost value (approx) 54.24\n')
(theta,J_history) = gradientDescent(X, y, theta, alpha, iterations)
print('Theta found by gradient descent:')
print(theta)
print('Expected theta values (approx):\n-3.6303\n 1.1664')
#plot data
plt.plot(X[:,1],np.dot(X,theta),c='b')
#predict
predict1 = np.dot([1,3.5],theta)
print('For population = 35,000, we predict a profit of %f' %(predict1 * 10000))
predict2 = np.dot([1,7],theta)
print('For population = 70,000, we predict a profit of %f' %(predict2 * 10000))
# Part 4: Visualizing J(theta_0, theta_1)
theta0_vals = np.linspace(-10,10,100)
theta1_vals = np.linspace(-1,4,100)
J_val = np.zeros((len(theta0_vals),len(theta1_vals)))
for i in range(1,len(theta0_vals)):
for j in range(1,len(theta1_vals)):
t = np.row_stack((theta0_vals[i],theta1_vals[j]))
J_val[i, j] = computeCost(X, y, t)
J_val = J_val.T
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(theta1_vals, theta0_vals, J_val,
rstride=1, cstride=1, cmap='rainbow')
plt.ylabel('theta_0')
plt.xlabel('theta_1')
plt.figure()
plt.contourf(theta0_vals, theta1_vals, J_val,
20, alpha=0.6, cmap=plt.cm.hot)
plt.plot(theta[0], theta[1], c='r', marker='x', markerSize=10, LineWidth=2)
运行结果：
Part1.
[[1. 0. 0. 0. 0.] [0. 1. 0. 0. 0.] [0. 0. 1. 0. 0.] [0. 0. 0. 1. 0.] [0. 0. 0. 0. 1.]] Part2.
y: (97, 1) Part3.
X: (97, 2) With theta = [0 ; 0],Cost computed = 32.072734 Expected cost value (approx) 32.07
With theta = [-1 ; 2],Cost computed = 54.242455 Expected cost value (approx) 54.24
Theta found by gradient descent: [[-3.62981201] [ 1.16631419]] Expected theta values (approx): -3.6303 1.1664 For population = 35,000, we predict a profit of 4522.876458 For population = 70,000, we predict a profit of 45343.872966
参考资料：
https://blog.csdn.net/lccflccf/category_8379707.html
https://blog.csdn.net/Cowry5/article/details/83302646
https://blog.csdn.net/weixin_44027820/category_9754493.html
总结：遇到bug或是计算结果相差太大，注意检查矩阵的维度，多打印出来看看
