数据建模与分析：线性回归小论文.doc

上海住房面积和房价的线性回归分析王明黔（上海大学 机电工程与自动化学院，上海200444） 摘要：在数据构建统计模型的学习中，统计学习是其一种基础的学习方法。本文针对城市人口数目与饮品连锁店利润的关系，就已有的数据进行线性回归分析，利用Matlab工具进行数据的线性回归模拟，进而得出城市人口数目与饮品连锁店利润的散点图、拟合直线图和三维等高线图。为了分析上海地区的住房面积和房价的关系，收集最近的售房成交数据，将数据导入到Matlab进行分析，得出上海房价与住房面积的线性关系。关键词：Matlab；线性回归；目标函数；梯度下降；统计学习基于数据的机器学习是现代智能技术中十分重要的一个方面，主要研究如何从一些观测数据(样本)出发，得出目前尚不能通过原理性分析得到的规律，并用以对未来数据或无法观测的数据进行预测。现实生活中大量存在我们尚无法准确认识但却可以进行观测的事件。因此，这种机器学习在从现代科学、技术到社会、经济等各领域中都有着十分重要的应用1。使用线性回归方法可以对一些观测数据进行分析，把预测事件中一些因素作为自变量，另一些随自变量变化而变化的变量作为因变量，研究它们之间的非确定性因果关系，以便预测因变量的未来发展趋势。根据若干观测数据寻找描述变量之间的函数或统计相关关系的最佳数学表达式，或者匹配数据之间相关关系的最佳拟合曲线，来表达随机性变量间的规律2。利用线性回归通过多变量机器学习的方法，可以建立上海住房面积和价格的线性关系，建立数学模型并评估其中的未知参数。1 案例分析1.1 目标函数的建立根据已知给出的城市人口数目与饮品连锁店利润的一些数据，可以得到一个样本集，如图1，为样本在Matlab软件加载数据图，第一列表示城市人口数目，第二列表示饮品连锁店利润。图1 城市人口数目与饮品连锁店利润的样本集Fig 1 Urban population and beverage chain profits of sample set针对样本集，我们可以假设一个线性模型： （1）式中：假设的线性模型；样本/特征；参数。其中，为模型参数，因此问题就可以转换为，求出的值。为了得到较为准确的回归线，应该满足一个前提，即各样本点尽可能分布在所建立的线性模型周围，因此我们建立目标函数： (2)式中：向量中的第个元素；向量中的第个元素；模型假设；训练集的数量。因此，我们只要求得使上述函数的值最小时的的值。1.2 目标函数的求解我们首先从一组开始，利用最速下降法不断改变的值来减小，直到达到我们希望得到的最小值3。对于最速下降法，可利用以下公式而来求解： (3)式中：学习速率（步长）； 下降方向。其中为修正学习速率，即搜索的步长，当过小，会延长搜索的时间，当过大，可能会错过极值点，导致不收敛。还有要注意最速下降法会收敛到局部最小值的情况，这种情况即便修正学习速率也不会改变结果，并且接近局部最小值时，因为斜率的减小最速下降法会自动减少每步的补偿4。因此，我们要现注意目标函数的极小值情况，如果有局部最小值可通过调节参数的大小来跳过局部最小值情况，否则，就不需跳过，最终通过多次迭代可以得到最终的理想结果。2 Matlab求解2.1 最速下降法的循环语句在Matlab程序中，最速下降法是一个不断迭代循环的过程，且确保同时更新，具体循环如下：2.2 ComputeCost函数的定义J函数即目标函数，m代表数据的个数，具体程序如下：function J = computeCost(X, y, theta)m = length(y);J = 0;predictions=X*theta;J=1/(2*m)*(predictions-y)*(predictions-y);end2.3 gradientDescent函数的定义GradientDescent函数就是最速下降法的迭代与循环过程5，具体如下：function theta, J_history = gradientDescent(X, y, theta, alpha, num_iters)m = length(y);J_history = zeros(num_iters, 1);for iter = 1:num_iterstemp1 = theta(1) - (alpha / m) * sum(X * theta - y).* X(:,1);temp2 = theta(2) - (alpha / m) * sum(X * theta - y).* X(:,2);theta(1) = temp1;theta(2) = temp2;J_history(iter) = computeCost(X, y, theta);endend2.4 散点图的绘制具体的数据已经根据一定格式记录在txt文件中，因此只需调取其中的数据并将散点绘制到图中6，具体程序如下：function plotData(x, y)figure; data = load(ex1data1.txt);X = data( : , 1 );Y = data( : , 2);X = ones(size(X,1),1),X; plot(X,Y,rx,MarkerSize, 4); axis(4 24 -5 25); xlabel(x); ylabel(y); end运行程序后得到散点图如图2所示：图2 散点图Fig.2 scatter diagram2.5 回归线的绘制通过机器统计学习后得到线性回归线，如图3所示： 图3 回归线fig.3 The regression line运行主程序中的语句计算出和J的最小值，程序如下：theta,J_history = gradientDescent(X, y, theta, alpha, num_iters);Matlab程序运行后结果报告如图4所示：图4 报告窗口fig.4 The report window其中ans结果为Matlab运行的最后一步结果，下面为模型参数值。假设模型参数值： ，2.6 等高线及三维图的绘制将和J的值绘制在三维图上，将的值以J为高度绘制绘制成等高线图，并将中心点即优化的参数用红色点标记，如下图5所示： 图5 三维图和等高线图Fig. 5 Graphic model and contour map等高线图中点红点代表了目标函数最小是对应的模型参数值，该参数值下的直线到每一个统计数据点的距离最短，最能代表离散统计数据的线性回归特性。2.7 结论根据已经给定的城市人口数目与饮品连锁店利润的关系，通过建立目标函数，由Matlab软件进行迭代求解，求解后的线性回归模型函数如下： (4)3 实际问题中国人多地少，土地不可再生，土地的稀缺性导致土地价格在未来城市发展中看涨，因此，房价也自然总体趋涨7。对于一些在沪求学的异地研究生们，在毕业后会存在是否留沪的选择。而房价是我们选择留沪的主要因素。为了了解上海的住房面积和价格的情况，我从链家网随机查找了120组不同区域以及不同面积的售房交易数据，将这些数据进行整理，并且希望通过统计学习的手段得出一般住房面积和价格的线性回归关系，通过设定两个期望住房面积来估计对应面积下的购房价格。3.1 程序的设计首先，通过EXCEL表格将收集的数据进行整理，得出了上海地区售房面积和价格对应表格，120组数据具体可见于附录。将120组数据通过TXT文件进行编辑后导入到Matlab软件中，其部分数据在软件中的表现形式如图6所示：图6 售房面积和价格数据在Matlab中表现形式Fig. 6 The data expression form of house area and price on sale in Matlab将这些数据以txt的格式保存，文件名为house，以便在主程序加载数据调用。该最速下降统计学习程序主要分为调用数据、绘制离散图、梯度下降和拟合直线、可视化显示目标函数、目标函数值与迭代次数关系函数这五个部分。3.1.1 调用数据程序将house.txt文件中的数据加载到data矩阵中，将data矩阵中的第一列中所有值组成一个向量储存在X向量中，将data矩阵中的第二列中所有值组成一个向量储存在y向量中，然后在y向量中计算训练样本数量存储到变量m中。其相关程序段如下：data = load(house.txt); %调用数据X = data(:, 1); y = data(:, 2); %向量存储m = length(y); %计算训练样本数量3.1.2 绘制离散图程序 该过程主要包括主程序的调用以及调用函数的定义，定义函数首先要打开一个绘图窗口，然后定义数值范围，最终附上标签绘制图形或者点。通过这样的方式处理，我们可以清晰的看到上海售房面积和价格的离散分布关系。主程序部分：plotData(X, y);fprintf(Program paused.Press enter to continue.n);定义调用函数plotData(X, y)function plotData(x, y)figure; % open a new figure windowplot(x, y, rx, MarkerSize,10);ylabel(Shanghai Commercial Housing Price in RMB one million);xlabel(Shanghai Commercial Housing Area in 10 square meters);end3.1.3 梯度下降和拟合直线程序此处先对参数进行初始化设置，然后设置初始的迭代次数和计算步长（学习速率），接着进行目标函数的迭代计算，并且不断将模型参数显示在屏幕上，接着画出拟合的直线，在屏幕上显示显示“数据训练和直线拟合”，最后根据拟合的直线估计80和160的商品房所对应的价格。主程序部分：fprintf(Running Gradient Descent .n)X = ones(m, 1), data(:,1); % Add a column of ones to xtheta = zeros(2, 1); % initialize fitting parameters iterations = 1000; % Some gradient descent settingsalpha =0.01; % compute and display initial costcomputeCost(X, y, theta) % run gradient descenttheta = gradientDescent(X, y, theta, alpha, iterations); % print theta to screenfprintf(Theta found by gradient descent: );fprintf(%f %f n, theta(1), theta(2); % Plot the linear fithold on; % keep previous plot visibleplot(X(:,2), X*theta, -)legend(Training data, Linear regression)hold off % dont overlay any more plots on this figurepredict1 = 1,8 *theta; % Predict values for house area of 35,000 and 70,000fprintf(For house area = 80, we predict a price of %fn,. predict1*1000000);predict2 = 1, 16 * theta;fprintf(For house area = 160, we predict a price of %fn,. predict2*1000000);fprintf(Program paused. Press enter to continue.n);pause;定义梯度下降调用函数：function theta, J_history = gradientDescent(X, y, theta, alpha, num_iters)m = length(y); % number of training examplesJ_history = zeros(num_iters, 1);for iter = 1:num_iterstemp1 = theta(1) - (alpha / m) * sum(X * theta - y).* X(:,1);%最速下降发的迭代循环 temp2 = theta(2) - (alpha / m) * sum(X * theta - y).* X(:,2); theta(1) = temp1; theta(2) = temp2;%同时更新theta1和theta2 J_history(iter) = computeCost(X, y, theta);endend3.1.4 可视化显示目标函数 该处Matlab程序主要将目标函数和两模型参数的对应关系通过空间平面和等高线的形式表达出来，根据空间平面我们可以清晰看出目标函数的具体极值分布形式。根据等高线图可以看出最终迭代计算结果是否为最优解，倘若红色的点位于等高线的最中间，则该情况为最优解，否则，不是最优解。主程序部分：fprintf(Visualizing J(theta_0, theta_1) .n)theta0_vals = linspace(-10, 10, 100);theta1_vals = linspace(-1, 4, 100);J_vals = zeros(length(theta0_vals), length(theta1_vals);for i = 1:length(theta0_vals) for j = 1:length(theta1_vals) t = theta0_vals(i); theta1_vals(j); J_vals(i,j) = computeCost(X, y, t); endendJ_vals = J_vals;figure;surf(theta0_vals, theta1_vals, J_vals)xlabel(theta_0); ylabel(theta_1);figure;contour(theta0_vals, theta1_vals, J_vals, logspace(-2, 3, 20)xlabel(theta_0); ylabel(theta_1);hold on;plot(theta(1), theta(2), rx, MarkerSize, 10, LineWidth, 10);定义目标函数： function J = computeCost(X, y, theta)m = length(y); J = 0;predictions=X*theta;J=1/(2*m)*(predictions-y)*(predictions-y);end3.2 确定计算步长我们在这里迭代计算时，并不是一开始就能迭代成功，正如前面所说步长过小，会延长搜索的时间，步长过大，可能会错过极值点，导致不收敛。故我们可通过建立目标函数值和迭代次数的关系函数，将不同的步长以一定的数量级排列，分别作出目标函数值和已运算迭代次数的关系函数曲线。为了更好的展示目标函数值和已运算迭代次数的对应关系，可以设置足够的迭代次数，设定足够的迭代次数为3000。目标函数值与迭代次数关系函数程序段：iterations=3000;alpha=1;%输入第一个步长theta,J1 = gradientDescent(X, y, zeros(2,1), alpha, iterations);alpha=2;%输入第二个步长theta,J2 = gradientDescent(X, y, zeros(2,1), alpha, iterations);alpha=3;%输入第三个步长theta,J3 = gradientDescent(X, y, zeros(2,1), alpha, iterations);figure;plot(1:numel(J1), J1, -k, LineWidth, 2);%第一条线显示连续线hold on;plot(1:numel(J2), J2, - -k, LineWidth, 2);%第二条线显示虚线hold on;plot(1:numel(J3), J3, -.k, LineWidth, 2);%第三条线显示点画线axis(0 3000 0 15);%限定坐标值显示范围xlabel(迭代次数); ylabel(J的值);legend(alpha=1 , alpha=2, alpha=3)hold off现在选取步长为1、0.1、0.01和0.001、0.0001、0.00001两组数据进行迭代运算，通过Matlab作出目标函数值和已运算迭代次数的关系函数曲线，判断目标函数值的收敛性。 3.2.1 第一组迭代运算判断敛散性将第一组步长值1、0.1和0.01带入目标函数值与迭代次数关系函数程序段中进行运算作图，为了让图形表示清晰，其坐标值可再一次限定其范围，其运行结果如图7所示。根据图中结果可知：步长为1和0.1时，J的函数是发散的，而步长为0.01时，J的函数是收敛的。图7 第一组步长中目标函数值的敛散性Fig.7 The convergence curves of first group alpha 3.2.2 第二组迭代运算判断敛散性 将第一组步长值0.001、0.0001和0.00001带入目标函数值与迭代次数关系函数程序段中进行运算作图，为了让图形表示清晰，其坐标值可再一次限定其范围，其运行结果如图8所示。根据图中结果可知：在该组所有步长下，J的函数都是收敛的。图8 第二组步长中目标函数值的敛散性Fig.8 The convergence curves of second group alpha3.2.3 综合步长分析在第一组步长分析过程中，当步长为0.01时，J的函数虽然判断是收敛的，但是其具体在多少次迭代下已经收敛到最小值我们不能从图7中看出。在第二组步长分析过程中，虽然3个步长下均能保证目标函数值能收敛，但是当步长为0.00001时，其所需的迭代次数需达到5000才收敛，其计算时间太长，故不能算是最速下降。为确定合适的步长，现将步长为0.01、0.001和0.0001进行目标函数收敛分析。如图9所示，该图在不同的迭代次数尺度下表现不同步长对应目标函数J的收敛特性。图9 综合目标函数收敛特性分析Fig.9 The convergence curves of composite alpha根据上图和Matlab迭代计算结果的最小值和对应最小值所需的最少迭代次数填入下表，如下表1。表1 步长和目标函数值Tab.1 The alpha and J步长alpha目标函数最小值达到最小值时所需最少迭代次数0.012.49592914980.0012.495929149880.00012.495929149892由上表可知选用步长数量级为0.01，可实现最速梯度下降。此处，不再讨论0.01附近的步长所对应的目标函数迭代次数关系。3.3 计算假设模型参数上面已经得到比较适合的迭代参数，现在可将步长为0.01，迭代次数iterations为1500带入到程序中进行计算，其运算结果如下：ans =14.2726Theta found by gradient descent: -0.715237 0.590743运算结果中ans代表了Matlab软件最后一步运算的结果，gradient descent代表了模型参数值，即模型参数值，即参数、。以下为可视化的点图、拟合直线图、目标函数与模型参数的立体图以及等高线图。通过plot函数作图可以得出售房面积和售房价格的离散数据图。如图10所示，从离散图可以看出数据整体是呈现售房价格与售房面积线性递增的趋势。图10 离散点图Fig. 10 The discrete point diagram接下来通过统计学习的方法拟合直线，如图11所示，从拟合直线可以看出拟合后的售房房价与售房面积的线性关系，直线以上数据差不多为市中心售房价位，直线以下大概为郊区售房价位，直线上大概为市中心与郊区的普遍价位。图11. 线性拟合图Fig11. Linear regression figure接下来可以通过surf函数和contour函数绘制目标函数与模型参数的立体图以及等高线图。如图12所示，从曲面图可以看出目标函数随模型参数的变化关系，可以明显看出该目标函数只有一个极小值，因此该处不会出现局部最小值情况。根据目标函数等高线图可以看出迭代计算后的参数值位于等高线的中间区域，因此可以判断出以上合适步长优化分析过程是正确的。图12 实际问题三维图和等高线图Fig. 12 Graphic model and contour map in practical problem通过上面的分析，这里的迭代过程显然是合理的，我们可以直接得出模型的线性回归方程：在程序中我们添加了预测情况计算，其具体主程序段如下：predict1 = 1,8 *theta;fprintf(For house area = 80, we predict a house price of %fn,. predict1*1000000);predict2 = 1, 120 * theta;fprintf(For house area = 120, we predict a house price of %fn,. predict2*1000000);程序运算结果解释如下，当我们考虑购买上海商品房面积为80平方和120平方时，根据所求的线性方程计算所得的房价为4010703.397050元和70173870.477535元。Matlab中具体输出结果如下：For house area = 80, we predict a house price of 4010703.397050For house area = 120, we predict a house price of 70173870.4775353.4 步长与迭代次数的可视化分析上面已经找到合适的步长和迭代次数，现在我们使用控制变量法对步长和迭代次数进一步分析，通过等高线图的可视化手段来描述步长与迭代次数的相互关系。3.4.1 相同步长不同迭代次数可视化分析 我们首先控制合适的步长为0.01，将迭代次数设置在1498次范围附近，分析迭代次数对最终迭代计算极小值点的位置。为了使等高线表达得更清晰，这里我们通过多次迭代计算得出相同幅度的等高线图，通过得出的模型参数的极值点的位置可以发现不同迭代次数下的极值点的移动规律。如图13所示，该图表现出四种不同迭代次数下的极点位置，根据这些位置，我们可以总结出极点位置特征。其特征总结如下表2所示。表2 迭代次数取值与梯度下降关系Tab.2 The relationship betweem iterations and gradient descent步长alpha迭代次数最终梯度下降情况0.011位于第二环与第三环之间0.0110紧靠第一环内侧0.01500第一环内侧靠中0.011500中心位置图13 相同步长下的等高线图Fig 13 contour map in the same alpha由上述可知，当步长一定，当迭代次数大于等于最小迭代次数时，其模型参数值的坐标点位于正中心，然后随着迭代次数的减少，坐标点会由内往外变动。3.4.2 相同迭代次数和不同步长的可视化分析我们首先选择合适的迭代次数为1500次，将步长（学习速率）设置在0.01次范围附近，分析不同步长对最终迭代计算极小值点的位置。为了使等高线表达得更清晰，这里我们通过多次迭代计算得出相同幅度的等高线图，通过得出的模型参数的极值点的位置可以发现不同步长下的极值点的移动规律。如图14所示，该图表现出六种不同步长下的极点位置，我们可以总结出极点位置特征。其特征总结如下表3所示。表3 迭代次数取值与梯度下降关系Tab.3 The relationship betweem iterations and gradient descent迭代次数步长alpha最终梯度下降情况15000.00001位于第六环上0.00002位于第二环与第三环之间0.0001紧靠第一环内侧0.001第一环内侧靠中0.01中心位置0.23中心位置图14 相同迭代次数下的等高线图Fig 14 contour map in the same iterations由上述可知，当迭代次数一定，当步长大于0.01且不超过0.023时，其模型参数值的坐标点位于正中心，然后随着迭代步长的减少，坐标点会由内往外变动。4 总结以上对两个案例进行了详细的分析，现将上述内容模型结果进行总结：1.城市人口数量和饮料连锁店利润的线性回归方程为：2.上海住房面积和售房价格的线性回归方程：通过统计线性回归的理论学习和Matlab实践，我充分了解了线性回归模型建立与求解的整个过程，学习了Matlab工具中最速下降法的迭代求解过程，熟悉使用Matlab软件的图形可视化的表示功能，为以后的课题数据的统计研究打下了良好的基础。致谢 感谢姚教授对本工作的大力支持，在此表示感谢！参考文献：1 张向君.信号分析与数据统计学习M.哈尔滨:哈尔滨工程大学出版社,2009:129-131.2 许继平,刘载文.同步变压位机械通气系统研究M.北京:北京大学出版社,2012:69-71.3 王玉英.优化与决策M.西安:西安交通大学出版社,2014:80-81.4 王桥.数字图像处理M.北京:科学出版社,2009:333-337.5 王正林,龚纯,何倩.精通Matlab科学计算M.北京:电子工业出版社,2012:401-404.6 刘浩.MatlabR2014a完全自学一本通M.北京:电子工业出版社,2015:1-708.7 郑友林.身边的经济学M.北京:中国市场出版社,2006:105-107.附录：表4 上海地区售房面积和价格表Tab.4 The house area and price on sale in Shanghai序号区域面积/价格/万元序号区域面积/价格/万元1黄浦 世博滨江51.5245544杨浦 中原75.333402松江大学城62.8524545闸北 阳城45.942653徐汇 上海南站52.424046松江 九亭134.884854松江 新桥51.7132247浦东 川沙85.962405虹口 江湾镇82.9245548浦东 唐镇87.663556杨浦 中原35.9515049浦东 航头93.542307嘉定 江桥39.495250青浦 盈浦63.82858浦东 金杨52.9626051松江 松江新城105.712589闵行 莘庄63.7139552闵行 吴泾21.795010长宁 北新泾36.7720053浦东 书院镇162.1413011松江 松江老城156.0131454浦东 外高桥74.7120812奉贤 西渡73.1815555浦东 川沙91.5228013普陀 长征56.9329056杨浦 中原69.8131014普陀 长征83.7141557闵行 古美76.8637515浦东 世博44.624558闵行 静安新城65.6634916宝山 淞南61.6123559闵行 老闵行128.5248017闵行 梅陇127.2258060闸北 西藏北路37.115918嘉定 江桥112.6740261浦东 祝桥60.420219浦东 金杨74.3439562普陀 光新136.6190020宝山 共富51.3716663闵行 古美121.3972021普陀 长寿路39.679064宝山 上大60.139022长宁 新华路59.5941065普陀 长征44.3422023浦东 三林103.0963066宝山 张庙8