多特征值的线性回归 梯度下降法

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

filepath='/Users/husir/Desktop/ex1data2.txt'
datafile=open(filepath, 'r')
data=datafile.readlines()
#print(*data)
#print(len(data))

size=[float(x.split(',')[0]) for x in data]
bednum=[float(x.split(',')[1]) for x in data]
price=[float(x.split(',')[2]) for x in data]

# 初步画个三维散点图
#ax=plt.subplot(111,projection='3d') # 创建一个三维的绘图工程

fig=plt.figure()
ax=Axes3D(fig)
ax.scatter(size,bednum,price)
ax.set_xlabel('size of house(in square feet)')
ax.set_ylabel('number of bedrooms')
ax.set_zlabel('price')

1.png

可以发现两组特征值相差较大，一个是几千的一个是个位数，因此需先进行特征缩放，以免梯度下降速率过慢

# 进行特征压缩再画个散点图
size=[x/1000 for x in size]
#print(size)
price=[x/100000 for x in price]

fig=plt.figure()
ax=Axes3D(fig)
ax.scatter(size,bednum,price)
ax.set_xlabel('size of house(in square feet)/1000')
ax.set_ylabel('number of bedrooms')
ax.set_zlabel('price/100000')

2.png

计算代价函数并用梯度下降法更新theta

x=np.ones((47,3)) # 47是样本数
y=np.ones((47,1))
for i in range(47):
    x[i][1]=size[i]
    x[i][2]=bednum[i]
    y[i][0]=price[i]
theta=[[1],[1],[1]]
iterations=1500 # 迭代次数
alpha=0.01 # 学习速率

'''
print(x)
print('---------------')
print(y)
print('---------------')
print(theta)
'''

for i in range(iterations):
    t=np.dot(x,theta)-y
    der0=sum(t)/47 # 这里不点乘x[:0]是因为都是1所以偷懒省掉了
    der1=sum(np.multiply(t,x[:,1].reshape(x[:,1].shape[0],1)))/47 # multiply是矩阵点乘
    der2=sum(np.multiply(t,x[:,1].reshape(x[:,1].shape[0],1)))/47
    theta[0]-=alpha*der0
    theta[1]-=alpha*der1
    theta[2]-=alpha*der2

print(*theta)
[-0.52910798] [0.76959084] [0.76959084]

作效果图

# 绘图
%matplotlib notebook # 为了能在jupyter notebook中3D旋转
fig=plt.figure()
ax=Axes3D(fig)

x=np.arange(0,5,0.1)
y=np.arange(1,5,1)
X,Y=np.meshgrid(x,y)
z=theta[0]+theta[1]*X+theta[2]*Y
ax.plot_surface(X,Y,z,cmap='coolwarm')

ax.scatter(size,bednum,price)
ax.set_xlabel('size of house(in square feet)/1000')
ax.set_ylabel('number of bedrooms')
ax.set_zlabel('price/100000')

3.png

补充

1. 矩阵点乘：np.multiply(a,b)
   矩阵乘法：np.dot(a,b) 或 np.matmul(a,b) 或 a.dot(b) 或 a@b

2. 提取矩阵的某行：mat[1]
   提取矩阵的某列：mat[:1]
   【注意！！！】提取矩阵某列的结果不会以n✖️1的矩阵返回，而会以1✖️n的返回。如这个例子中x[:,1]如果不reshape，就会得到这样的结果：

   
   
   4.png

x[:,1].reshape(x[:,1].shape[0],1)的作用就是将该矩阵重新塑造为n*1的向量。

3. 矩阵的转置：
   对于大于一维的矩阵：a.T 或 a.transpose()
   对于一维矩阵（向量）：a.reshape(a.shape[0],1)

上午都没搞物理唉 下午加油⛽️