本节是python实现多元回归的代码部分，理论参考链接: link.
代码下载地址link.
代码可直接赋值运行，如有问题请留言
本节使用的数据是收入与年龄，性别关系的多元线性回归
1 基本环境设置
import numpy as np
import matplotlib.pyplot as pl
import matplotlib
matplotlib.rcParams['font.sans-serif']='SimHei' #画图正常显示中文
matplotlib.rcParams['font.family']='sans-serif'
matplotlib.rcParams['axes.unicode_minus']=False2 读取数据加载
def loadDataset(filename):
X=[]
Y=[]
with open(filename,'rb') as f:
for idx,line in enumerate(f):
line=line.decode('utf-8').strip()
if not line:
continue
eles=line.split(',')
if idx==0:
numFea=len(eles)
eles=list(map(float,eles))#map返回一个迭代对象
X.append(eles[:-1])
Y.append([eles[-1]])
return np.array(X),np.array(Y)预览下数据，数据是如下图所示
第一列性别，第二列年龄，第三列收入
3 y预估方法与误差方法设计
y 估计方法
def h(theta,X):
return np.dot(X,theta)误差和估计
def J(theta,X,Y):
return np.sum(np.dot((h(theta,X)-Y).T,(h(theta,X)-Y))/(2*m))
4 梯度下降设计
def bgd(alpha,maxloop,epsilon,X,Y):
m,n=X.shape
theta=np.zeros((n,1))
count=0
converged=False
error=np.inf
errors=[J(theta,X,Y),]
thetas={}
for i in range(n):
thetas[i]=[theta[i,0],]
while count<=maxloop:
if(converged):
break
count=count+1
for j in range(n):
deriv=np.dot(X[:,j].T,(h(theta,X)-Y)).sum()/m
thetas[j].append(theta[j,0]-alpha*deriv)
for j in range(n):
theta[j,0]=thetas[j][-1]
error=J(theta,X,Y)
errors.append(error)
if(abs(errors[-1]-errors[-2])<epsilon):
converged=True
return theta,errors,thetas
5 数据处理
这里的数据没有异常值，缺失值。在R部分也讲过缺失值核异常值的处理，盖帽法填补，删除，或spss回归，knn填补
def standarize(X):
"""特征标准化处理
Args:
X 样本集
Returns:
标准化后的样本集
"""
m,n=X.shape
#归一化每一个特征
for j in range(n):
features=X[:,j]
meanVal=features.mean(axis=0)
std=features.std(axis=0)
if std!=0:
X[:,j]=(features-meanVal)/std
else:
X[:,j]=0
return X
读取属于与预览维度
ori_X,Y=loadDataset('./data/income.csv')
print(ori_X.shape)
print(Y.shape)
结果如下：
(100, 2)
(100, 1)
6 模型运行
m,n=ori_X.shape
X=standarize(ori_X.copy())
X=np.concatenate((np.ones((m,1)),X),axis=1)
alpha=0.3
maxloop=5000
epsilon=0.0000000000000001
result=bgd(alpha,maxloop,epsilon,X,Y)
theta,errors,thetas=result
print(errors)
print(theta)
结果如下：
[24.33730066195, 13.505993103864864, 8.227069671967811, 5.646107653055884, 4.380348621596483, 3.7577018007906924, 3.450468610785003, 3.2983894621721754, 3.2228610964245137, 3.1852188396658816, 3.166387845485983, 3.1569292252947605, 3.152157461763612, 3.149738790112871, 3.1485065999017103, 3.1478754430280538, 3.1475502771380808, 3.147381731762091, 3.147293811506613, 3.1472476466986463, 3.14722324380204, 3.1472102569749523, 3.147203299011644, 3.1471995464497224, 3.1471975096191853, 3.1471963972545858, 3.1471957862237945, 3.147195448748745, 3.1471952614181884, 3.147195156950363, 3.147195098447062, 3.147195065560232, 3.1471950470106522, 3.1471950365163703, 3.1471950305635144, 3.147195027178879, 3.1471950252505327, 3.1471950241499296, 3.147195023520792, 3.147195023160675, 3.1471950229543095, 3.147195022835933, 3.147195022767972, 3.1471950227289267, 3.147195022706478, 3.1471950226935674, 3.1471950226861374, 3.147195022681859, 3.1471950226793957, 3.1471950226779772, 3.1471950226771592, 3.147195022676689, 3.147195022676417, 3.14719502267626, 3.1471950226761694, 3.1471950226761174, 3.147195022676087, 3.1471950226760703, 3.147195022676061, 3.1471950226760543, 3.147195022676051, 3.1471950226760494, 3.147195022676048, 3.147195022676048]
[[ 6.142094 ]
[ 2.16407412]
[-0.03431546]]
7 模型可视化
使用的是三维绘图
%matplotlib
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
import matplotlib.ticker as mtick
fittingFig=pl.figure(figsize=(16,12))
title='bgd:rate=%.3f,maxloop=%d,epsilon=%.3f \n'%(alpha,maxloop,epsilon)
ax=fittingFig.gca(projection='3d')
xx=np.linspace(0,1,100)
yy=np.linspace(0,100,100)
zz=np.zeros((100,100))
for i in range(100):
for j in range(100):
normalizegender=(xx[i]-ori_X[:,0].mean(0))/ori_X[:,0].std(0)
normalizeAge=(yy[j]-ori_X[:,1].mean(0))/ori_X[:,1].std(0)
x=np.matrix([[1,normalizegender,normalizeAge]])
zz[i,j]=h(theta,x)
xx,yy=np.meshgrid(xx,yy)
ax.zaxis.set_major_formatter(mtick.FormatStrFormatter('%.2e'))
ax.plot_surface(xx,yy,zz,rstride=1,cstride=1,cmap=cm.rainbow,alpha=0.1,antialiased=True)
xs=ori_X[:,0].flatten()
ys=ori_X[:,1].flatten()
zs=Y[:,0].flatten()
ax.scatter(xs,ys,zs,c='b',marker='o')
ax.set_xlabel(u'性别')
ax.set_ylabel(u'年龄')
ax.set_zlabel(u'收入')可以发现模型的平面将数据按照维度较好拟合
8 误差绘图
%matplotlib inline
errorsFig=pl.figure()
ax=errorsFig.add_subplot(111)
ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.2e'))
pl.plot(range(len(errors)),errors)
pl.xlabel(u'迭代次数')
pl.ylabel(u'代价函数')
pl.show()
链接: [link](https://download.csdn.net/download/weixin_44498127/85285750
).
