梯度下降和线性回归（二）波士顿房价python代码实现

一.多元线性回归方程

假设样本中有m个特征量，那么对应的线性回归方程如下
$\hat{y} =\theta_{1} x_{1} + \theta_{2} x_{2} +...+ \theta_{m} x_{m}$

二.损失函数的构造

假设样本中有n个训练集

训练集

Loss=\frac{1}{2n} \sum_{i=1}^n(y_{i} -\hat{y}_ {i} )

在m个输出变量y中，用实际值减去回归方程中的预测值，平方求和再求均值来反映回归方程中的输出值与实际输出值的偏差程度平方求和可以避免

y_{i} -\hat{y}_ {i}

不同的正负情况在求和过程中抵消偏差，假如对于某一列中

y_{i} -\hat{y}_ {i}

,另一列中

y_{i} -\hat{y}_ {i}

,实际这两列都有偏差，不平方的的话对应损失为0，显然不妥。

三.迭代过程

利用梯度下降原理来将损失函数不断变小直至收敛，如果对梯度下降不是很清楚的话，可以了解梯度下降和线性回归（一）附python代码实现 - 简书。
$\theta _{i} =\theta _{i} -\alpha \frac{\sigma L(\theta )}{\sigma \theta_{i} }$
$b=b-\alpha \frac{\sigma L(\theta )}{\sigma b }$
$\hat{y}_{i} =\theta_{1} x_{i1} + \theta_{2} x_{i2} +...+ \theta_{m} x_{im}$
$Loss=\frac{1}{2n} \sum_{i=1}^n(y_{i} -\hat{y}_ {i} )$
$\frac{\sigma L(\theta )}{\sigma \theta_{i} } =-\frac{1}{n} \sum_{i=1}^n(y_{i} -\hat{y}_{i} )\frac{1}{n} \sum_{i=1}^n x_{i}$
$\frac{\sigma L(\theta )}{\sigma b } =-\frac{1}{n} \sum_{i=1}^n(y_{i} -\hat{y}_{i} )$
所以：
$\theta _{i} =\theta _{i} +\frac{1}{n} \sum_{i=1}^n(y_{i} -\hat{y}_{i} )\frac{1}{n} \sum_{i=1}^n x_{i}$
$b=b+\alpha \frac{1}{n} \sum_{i=1}^n(y_{i} -\hat{y}_{i} )$

四.利用多元线性回归分析波士顿房价

#调用库

import matplotlib.pyplot as plt

import numpy as np

import pandas as pd

import math

#查看数据集的前30行

df = pd.read_csv("波士顿房价2.csv", index_col=None)

df.head(30)  # 查看数据集的前30行

效果如下：

训练集

需要下载数据集的可以访问百度网盘请输入提取码,提取码1p23

CRIM=floor.loc[:,'CRIM'].values
ZN=floor.loc[:,'ZN'].values
INDUS=floor.loc[:,'INDUS'].values
CHAS=floor.loc[:,'CHAS'].values
NOX=floor.loc[:,'NOX'].values
RM=floor.loc[:,'RM'].values
AGE=floor.loc[:,'AGE'].values
DIS=floor.loc[:,'DIS'].values
RAD=floor.loc[:,'RAD'].values
TAX=floor.loc[:,'TAX'].values
PTRATIO=floor.loc[:,'PTRATIO'].values
LSTAT=floor.loc[:,'LSTAT'].values
MEDV=floor.loc[:,'MEDV'].values

各个特征量含义如下：
CRIM: 城镇人均犯罪率
ZN: 住宅用地所占比例
INDUS: 城镇中非住宅用地所占比例
CHAS: 虚拟变量,用于回归分析
NOX: 环保指数
RM: 每栋住宅的房间数
AGE: 1940 年以前建成的自住单位的比例
DIS: 距离 5 个波士顿的就业中心的加权距离
RAD: 距离高速公路的便利指数
TAX: 每一万美元的不动产税率
PTRATIO: 城镇中的教师学生比例
B: 城镇中的黑人比例
LSTAT: 地区中有多少房东属于低收入人群
MEDV: 自住房屋房价中位数（也就是均价）
在进行梯度下降前，我们需要分割一下数据集，查看数据集大小后，按8:2的比例，将数据集的前406行作为训练集，后面100行作为测试集

训练集大小

测试集划分

测试集

受能力影响，笔者用最比较繁琐的代码实现了梯度下降更新，如果熟悉矩阵乘法，需要简单的实现，可以参考多元线性回归-波士顿房价预测问题python_W_yu_cheng的博客-CSDN博客_波士顿房价问题数学建模。代码如下

#梯度下降法优化拟合方程
#设定参数
learning_rate = 0.001
l=len(MEDV)-106
#损失函数
def L(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y,theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12):
    return pow(theta1*x1+theta2*x2+theta3*x3+theta4*x4+theta5*x5+theta6*x6+theta7*x7+theta8*x8+theta9*x9+theta10*x10+theta11*x11+theta12*x12+theta0-y,2)/(2*l)
#损失函数关于theta1求偏导
def L1(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y,theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12):
    return -x1*(theta1*x1+theta2*x2+theta3*x3+theta4*x4+theta5*x5+theta6*x6+theta7*x7+theta8*x8+theta9*x9+theta10*x10+theta11*x11+theta12*x12+theta0-y)/l
#损失函数关于theta2求偏导
def L2(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y,theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12):
    return -x2*(theta1*x1+theta2*x2+theta3*x3+theta4*x4+theta5*x5+theta6*x6+theta7*x7+theta8*x8+theta9*x9+theta10*x10+theta11*x11+theta12*x12+theta0-y)/l
#损失函数关于theta3求偏导
def L3(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y,theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12):
    return -x3*(theta1*x1+theta2*x2+theta3*x3+theta4*x4+theta5*x5+theta6*x6+theta7*x7+theta8*x8+theta9*x9+theta10*x10+theta11*x11+theta12*x12+theta0-y)/l
#损失函数关于theta4求偏导
def L4(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y,theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12):
    return -x4*(theta1*x1+theta2*x2+theta3*x3+theta4*x4+theta5*x5+theta6*x6+theta7*x7+theta8*x8+theta9*x9+theta10*x10+theta11*x11+theta12*x12+theta0-y)/l
#损失函数关于theta5求偏导
def L5(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y,theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12):
    return -x5*(theta1*x1+theta2*x2+theta3*x3+theta4*x4+theta5*x5+theta6*x6+theta7*x7+theta8*x8+theta9*x9+theta10*x10+theta11*x11+theta12*x12+theta0-y)/l
#损失函数关于theta6求偏导
def L6(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y,theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12):
    return -x6*(theta1*x1+theta2*x2+theta3*x3+theta4*x4+theta5*x5+theta6*x6+theta7*x7+theta8*x8+theta9*x9+theta10*x10+theta11*x11+theta12*x12+theta0-y)/l
#损失函数关于theta7求偏导
def L7(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y,theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12):
    return -x7*(theta1*x1+theta2*x2+theta3*x3+theta4*x4+theta5*x5+theta6*x6+theta7*x7+theta8*x8+theta9*x9+theta10*x10+theta11*x11+theta12*x12+theta0-y)/l
#损失函数关于theta8求偏导
def L8(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y,theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12):
    return -x8*(theta1*x1+theta2*x2+theta3*x3+theta4*x4+theta5*x5+theta6*x6+theta7*x7+theta8*x8+theta9*x9+theta10*x10+theta11*x11+theta12*x12+theta0-y)/l
#损失函数关于theta9求偏导
def L9(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y,theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12):
    return -x9*(theta1*x1+theta2*x2+theta3*x3+theta4*x4+theta5*x5+theta6*x6+theta7*x7+theta8*x8+theta9*x9+theta10*x10+theta11*x11+theta12*x12+theta0-y)/l
#损失函数关于theta10求偏导
def L10(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y,theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12):
    return -x10*(theta1*x1+theta2*x2+theta3*x3+theta4*x4+theta5*x5+theta6*x6+theta7*x7+theta8*x8+theta9*x9+theta10*x10+theta11*x11+theta12*x12+theta0-y)/l
#损失函数关于theta11求偏导
def L11(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y,theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12):
    return -x11*(theta1*x1+theta2*x2+theta3*x3+theta4*x4+theta5*x5+theta6*x6+theta7*x7+theta8*x8+theta9*x9+theta10*x10+theta11*x11+theta12*x12+theta0-y)/l
#损失函数关于theta12求偏导
def L12(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y,theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12):
    return -x12*(theta1*x1+theta2*x2+theta3*x3+theta4*x4+theta5*x5+theta6*x6+theta7*x7+theta8*x8+theta9*x9+theta10*x10+theta11*x11+theta12*x12+theta0-y)/l
#损失函数关于theta0求偏导
def L13(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y,theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12):
    return -(theta1*x1+theta2*x2+theta3*x3+theta4*x4+theta5*x5+theta6*x6+theta7*x7+theta8*x8+theta9*x9+theta10*x10+theta11*x11+theta12*x12+theta0-y)/l
#损失函数值变化记录数组
Loss=[]
#未迭代时的损失函数
Loss0=0
for i in range(l):
    Loss0=Loss0+L(CRIM[i],ZN[i],INDUS[i],CHAS[i],NOX[i],RM[i],AGE[i],DIS[i],RAD[i],TAX[i],PTRATIO[i],LSTAT[i],MEDV[i],10,10,10,10,10,10,10,10,10,10,10,10,10)
Loss.append(Loss0)
#初始化theta0 theta1
theta0=10
theta1=10
theta2=10
theta3=10
theta4=10
theta5=10
theta6=10
theta7=10
theta8=10
theta9=10
theta10=10
theta11=10
theta12=10
#进行迭代
for i in range(2000):
    dertheta0=0
    dertheta1=0
    dertheta2=0
    dertheta3=0
    dertheta4=0
    dertheta5=0
    dertheta6=0
    dertheta7=0
    dertheta8=0
    dertheta9=0
    dertheta10=0
    dertheta11=0
    dertheta12=0
    dertheta13=0
    Loss1=0
    for j in range(l):
        dertheta0=dertheta0+L13(CRIM[j],ZN[j],INDUS[j],CHAS[j],NOX[j],RM[j],AGE[j],DIS[j],RAD[j],TAX[j],PTRATIO[j],LSTAT[j],MEDV[j],theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12)
        dertheta1=dertheta0+L1(CRIM[j],ZN[j],INDUS[j],CHAS[j],NOX[j],RM[j],AGE[j],DIS[j],RAD[j],TAX[j],PTRATIO[j],LSTAT[j],MEDV[j],theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12)
        dertheta2=dertheta0+L2(CRIM[j],ZN[j],INDUS[j],CHAS[j],NOX[j],RM[j],AGE[j],DIS[j],RAD[j],TAX[j],PTRATIO[j],LSTAT[j],MEDV[j],theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12)
        dertheta3=dertheta0+L3(CRIM[j],ZN[j],INDUS[j],CHAS[j],NOX[j],RM[j],AGE[j],DIS[j],RAD[j],TAX[j],PTRATIO[j],LSTAT[j],MEDV[j],theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12)
        dertheta4=dertheta0+L4(CRIM[j],ZN[j],INDUS[j],CHAS[j],NOX[j],RM[j],AGE[j],DIS[j],RAD[j],TAX[j],PTRATIO[j],LSTAT[j],MEDV[j],theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12)
        dertheta5=dertheta0+L5(CRIM[j],ZN[j],INDUS[j],CHAS[j],NOX[j],RM[j],AGE[j],DIS[j],RAD[j],TAX[j],PTRATIO[j],LSTAT[j],MEDV[j],theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12)   
        dertheta6=dertheta0+L6(CRIM[j],ZN[j],INDUS[j],CHAS[j],NOX[j],RM[j],AGE[j],DIS[j],RAD[j],TAX[j],PTRATIO[j],LSTAT[j],MEDV[j],theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12)
        dertheta7=dertheta0+L7(CRIM[j],ZN[j],INDUS[j],CHAS[j],NOX[j],RM[j],AGE[j],DIS[j],RAD[j],TAX[j],PTRATIO[j],LSTAT[j],MEDV[j],theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12)  
        dertheta8=dertheta0+L8(CRIM[j],ZN[j],INDUS[j],CHAS[j],NOX[j],RM[j],AGE[j],DIS[j],RAD[j],TAX[j],PTRATIO[j],LSTAT[j],MEDV[j],theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12) 
        dertheta9=dertheta0+L9(CRIM[j],ZN[j],INDUS[j],CHAS[j],NOX[j],RM[j],AGE[j],DIS[j],RAD[j],TAX[j],PTRATIO[j],LSTAT[j],MEDV[j],theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12)
        dertheta10=dertheta0+L10(CRIM[j],ZN[j],INDUS[j],CHAS[j],NOX[j],RM[j],AGE[j],DIS[j],RAD[j],TAX[j],PTRATIO[j],LSTAT[j],MEDV[j],theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12)
        dertheta11=dertheta0+L11(CRIM[j],ZN[j],INDUS[j],CHAS[j],NOX[j],RM[j],AGE[j],DIS[j],RAD[j],TAX[j],PTRATIO[j],LSTAT[j],MEDV[j],theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12)
        dertheta12=dertheta0+L12(CRIM[j],ZN[j],INDUS[j],CHAS[j],NOX[j],RM[j],AGE[j],DIS[j],RAD[j],TAX[j],PTRATIO[j],LSTAT[j],MEDV[j],theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12)
    theta0=theta0+learning_rate*dertheta0
    theta1=theta1+learning_rate*dertheta1
    theta2=theta2+learning_rate*dertheta2
    theta3=theta3+learning_rate*dertheta3
    theta4=theta4+learning_rate*dertheta4
    theta5=theta5+learning_rate*dertheta5
    theta6=theta6+learning_rate*dertheta6
    theta7=theta7+learning_rate*dertheta7
    theta8=theta8+learning_rate*dertheta8
    theta9=theta9+learning_rate*dertheta9
    theta10=theta10+learning_rate*dertheta10
    theta11=theta11+learning_rate*dertheta11
    theta12=theta12+learning_rate*dertheta12
    for k in range(l):
        Loss1=Loss1+L(CRIM[k],ZN[k],INDUS[k],CHAS[k],NOX[k],RM[k],AGE[k],DIS[k],RAD[k],TAX[k],PTRATIO[k],LSTAT[k],MEDV[k],theta0,theta1,theta2,theta3,theta4,theta5,theta6,theta7,theta8,theta9,theta10,theta11,theta12)
    Loss.append(Loss1)

在梯度下降更新完成后，我们来看一下损失函数的收敛曲线图

损失函数收敛图

将特征变量放入一个数组中

特征变量数组

最后我们来查看一下测试集情况

for i in range(100):

    print("true:\t{}".format(test_y[i]),end="\t")

    pre = np.dot(theta,test_x[i])+theta0

    print("guess:\t{}".format(pre))

plt.show()

测试集情况

最后编辑于：2022.01.27 20:47:31

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 217,734评论 6赞 505
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 92,931评论 3赞 394
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 164,133评论 0赞 354
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 58,532评论 1赞 293
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 67,585评论 6赞 392
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 51,462评论 1赞 302
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 40,262评论 3赞 418
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 39,153评论 0赞 276
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 45,587评论 1赞 314
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 37,792评论 3赞 336
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 39,919评论 1赞 348
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 35,635评论 5赞 345
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 41,237评论 3赞 329
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 31,855评论 0赞 22
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,983评论 1赞 269
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 48,048评论 3赞 370
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 44,864评论 2赞 354