还是校数学建模。。。发了N篇博客。。。这个代码是判断影响因素权重比的。 
 今天总算学会怎么把代码插入了。。。
 
n=10;m=4;
x1=[ 760.72 773.48 784.17 794.62 806.14 814.58 822.30 832.31 842.42 854.19];
x2=[ 36770 40561 45702 49518 54494 57474 67515 73678 74246 81171];
x3=[ 6.62 7.49 6.99 6.99 6.27 6.82 6.68 6.68 6.15 5.40];
x4=[7831 8831 9870 12754 14021 16120 16672 19753 22198 25407];
y=[6795,7455.6,8677.6,11534.9,14605.5,15472.1,18163.9,18757.1,19711.5,21163.5];

X=[ones(n,1), x1',x2',x3',x4'];
[b,bint,r,rint,s]=regress(y',X);
s2=sum(r.^2)/(n-m-1);
%b,bint,s,s2;
rcoplot(r,rint);

1. Matlab函数
 函数：regress [b,bint,r,rint,stats]=regress(y,x);
 参数：x、y
 X：矩阵x表示预测变量，要注意的是如果有常数项，要在矩阵x中包含一个由1构成的列。通常的处理是 x = [ones(n,1),x1] x1表示原有的预测变量，n表示数据的组数。
 返回值：b、bint、r、rint、stats
 b：多元线性回归的系数估计值
 bint：还返回系数估计值的 95% 置信区间的矩阵 bint。
 r ：还返回由残差组成的向量 r。
 rint ：返回矩阵 rint，其中包含可用于诊断离群值的区间。
 Stats:返回向量 stats，其中包含 R2 统计量、F 统计量及其 p 值，以及误差方差的估计值。
 矩阵 X 必须包含一个由 1 组成的列，以便软件正确计算模型统计量。
 
2.实现简单的一元线性回归
 例1. 合金的强度 y 与其中的碳含量 x 有比较密切的关系，今从生产中收集了一批数据如下表 1
 %x=[0.1:0.01:0.18]’;
 %y=[42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0]’;
 
 
x=0.1:0.01:0.18; y=[42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0]; x1=[0.1:0.01:0.18]'; y=[42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0]'; x=[ones(9,1),x1]; 
[b,bint,r,rint,stats]=regress(y,x); 
x=0.1:0.01:0.18; y=[42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0]; b=b';x_0=0.1:0.01:0.18;
 y_0=b(2)*x+b(1); 
 plot(x,y,'+',x_0,y_0)
 

 
 
分析：beta0=27.4722，beta1=137.5000
 置信区间是beta0[18.6851,36.2594]，beta1[75.7755,199.2245]；R2=0.7985,F = 27.7469, p = 0.0012,S2=4.0883。
 
3. 以此类推到多元线性回归
 例 2 某厂生产的一种电器的销售量 y 与竞争对手的价格 1 x 和本厂的价格 2 x 有关。表 2 是该商品在 10 个城市的销售记录。试根据这些数据建立 y 与 1 x 和 2 x 的关系式，对得到的模型和系数进行检验。若某市本厂产品售价 160（元），竞争对手售价 170（元），预测商品在该市的销售量。
 
x1=[120 140 190 130 155 175 125 145 180 150]';
 x2=[100 110 90 150 210 150 250 270 300 250]'; 
 y=[102 100 120 77 46 93 26 69 65 85]'; 	x=[ones(10,1),x1,x2];
  [b,bint,r,rint,stats]=regress(y,x); b,bint,stats

二、多元线性回归原理
 
2.1、数学模型
 
在社会生活及生产实践中会经常遇到一种问题，即我们非常关注一个量的变化，而这个量受到另一个或是多个因素的影响，我们想要了解这些因素是如何影响我们最为关注的这个量的以及这些因素对我们最为关注的这个量的影响权重分别有多大，知道了这些，我们就可以对该量变化所反映的相关问题做出分析和评价，并对其未来发展趋势进行预测和控制，这里就要用到数理统计中一个非常重要而普遍的分析方法，即回归分析法。
 
如果一个因变量y与k个自变量$x_1,x_2,\dots ,x_k$存在线性相关关系，那么就可以用多元线性回归模型
 $$y=a_0+a_1{x}_1+a_2{x}_2+\dots +a_k{x}_k,\text{公式1}$$
 对其进行描述，其中未知常量$a_0,a_1,\dots ,a_k$称为回归模型系数，若n次抽样，第$i$次抽样数据为$\left(y_i,x_{1i},x_{2i},\cdots ,x_{ki}\right)$那么就有
 $$\{\begin{array}{l}{y}_1=a_0+a_1{x}_{11}+a_2{x}_{21}+\dots +a_k{x}_{k1}+{\epsilon }_1\\ y_2=a_0+a_1{x}_{12}+a_2{x}_{22}+\dots +a_k{x}_{k2}+{\epsilon }_2\\ ⋮\\ y_n=a_0+a_1{x}_{1i}+a_2{x}_{2i}+\dots +a_i{x}_{ki}+{\epsilon }_n\end{array}\text{公式2}$$
 其中${\epsilon }_0,{\epsilon }_1,\dots ,{\epsilon }_n$为随机误差项，回归分析的主要任务就是以误差${\epsilon }_0,{\epsilon }_1,\dots ,{\epsilon }_n$的平方和最小为原则，求多元回归模型的回归系数$a_0,a_1,\dots ,a_k$。
 
求解这个方程是要以$S={\sum }_{i=1}^{i=n}{\epsilon }_i^2={\sum }_{i=1}^{i=n}{\left(a_0+a_1{x}_{1i}+\cdots +a_k{x}_{ki}-y_i\right)}^2$为最小原则，求$a_0,a_1,\dots ,a_k$要使得S最小，应该满足$\frac{\partial S}{\partial a_j}=0,j=0,1,\dots ,k$
 $$即\{\begin{array}{l}{\sum }_{i=1}^{i=n}2\left(a_0+a_1{x}_{1i}+a_2{x}_{2i}+\cdots +a_k{x}_{ki}-y_i\right)=0\\ {\sum }_{i=1}^{i=n}2\left(a_0+a_1{x}_{1i}+a_2{x}_{2i}+\cdots +a_k{x}_{ki}-y_i\right)x_{1i}=0\\ ⋮\\ {\sum }_{i=1}^{i=n}2\left(a_0+a_1{x}_{1i}+a_2{x}_{2i}+\cdots +a_k{x}_{ki}-y_i\right)x_{ni}=0\end{array},\text{公式3}$$
 
$$有\{\begin{array}{l}n{a}_0+{\sum }_{i=1}^{i=n}{x}_{1i}{a}_1+\cdots +{\sum }_{i=1}^{i=n}{x}_{ki}{a}_k={\sum }_{i=1}^{i=n}{y}_i\\ {\sum }_{i=1}^{i=n}{x}_{1i}{a}_0+{\sum }_{i=1}^{i=n}{x}_{1i}^2{a}_1+\cdots +{\sum }_{i=1}^{i=n}{x}_{11}{x}_{ki}{a}_k={\sum }_{i=1}^{i=n}{x}_{1i}{y}_i\\ ⋮\\ {\sum }_{i=1}^{i=n}{x}_{ki}{a}_0+{\sum }_{i=1}^{i=n}{x}_{1i}{x}_{ki}{a}_1+\cdots +{\sum }_{i=1}^{i=n}{x}_{ki}^2{a}_k={\sum }_{i=1}^{i=n}{x}_{ki}{y}_i\end{array},\text{公式4}$$
 
上式可以写成形式：$Y=XA$
 
其中：
 $$X=\left[\begin{array}{cccc}n& {\sum }_{i=1}^{i=n}{x}_{1i}& \dots & {\sum }_{i=1}^{i=n}{x}_{ki}\\ {\sum }_{i=1}^{i=n}{x}_{1i}& {\sum }_{i=1}^{i=n}{x}_{1i}^2& \dots & {\sum }_{i=1}^n{x}_{11}{x}_{ki}\\ ⋮& ⋮& \dots & ⋮\\ {\sum }_{i=1}^{i=n}{x}_{ki}& {\sum }_{i=1}^{i=n}{x}_{11}{x}_{ki}& \dots & {\sum }_{i=1}^{i=n}{x}_{ki}^2\end{array}\right],Y=\left[\begin{array}{c}{\sum }_{i=1}^{i=n}{y}_i\\ {\sum }_{i=1}^{i=n}{x}_{1i}{y}_i\\ ⋮\\ {\sum }_{i=1}^{i=n}{x}_{2i}{y}_i\end{array}\right],A=\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_k\end{array}\right]$$
 
公式2也可以直接写成矩阵表达式：$Y=XA+E$
 
其中：
 $$Y=\left[\begin{array}{l}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right],X=\left[\begin{array}{ccccc}1& x_{1,1}& \dots & x_{k,1}& \\ 1& x_{1,2}& \dots & x_{k,2}& \\ ⋮& ⋮& \dots & ⋮& \\ 1& x_{1,k}& \dots & x_{2,k}& \end{array}\right],A=\left[\begin{array}{l}{a}_0\\ a_1\\ ⋮\\ a_k\end{array}\right],E=\left[\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\\ ⋮\\ {\epsilon }_n\end{array}\right]$$
 那么我们的任务就是求解$A$
 
2.2、案例分析
 
分析总能耗与其他变量的关系，首先分析其相关性，做相似矩阵分析，如下：
 
 
第一步：得出总能耗与那些变量之间的关系，求出相关系数
 
 
[Data,str]=xlsread('C:\Users\86188\Desktop\仿真数据\北京参数一百组新.xlsx','sheet1','A1:R101',0.4);% 得到表格中所有数据
Resemblance=corrcoef(Data);                                 % 得到系数相关矩阵
[Row,Col]=size(Resemblance);                                % 得到Resemblance矩阵的行和列
site=[];count=0;                                            % 
for row=1:(Row-1)                                           % 得到总能源与哪些因素有关
    if(abs(Resemblance(row,Col))>0.4)
        disp(['总能源与第' num2str(row) '列' str(row) '相关度较高，相关系数为 ' num2str(Resemblance(Col,row)) ]);
        count=count+1;                                      % 保存相关变量的数量
        site(count)=row;                                    % 保存相关变量的列地址
    end
end
 
第二步：得到回归系数和置信区间
 
 
[Row,Col]=size(Data);                                       % 得到数据矩阵的行和列的大小
ConVariable=zeros(Row,count);                               % 创建一个行相等列指定的矩阵                 
for row=1:count
    ConVariable(:,row)=Data(:,site(row));
end
TotalEnergy=Data(:,Col);                                   % 得到总能源原始数据
IndeVariable=[ones(Row,1),ConVariable];                    % 创建相关变量数组
[b,bint,r,rint,stats]=regress(TotalEnergy,IndeVariable);   % 求回归系数的点估计和区间估计、并检验回归模型
 
第三步：做残差分析，得出那些点偏差太大
 
 
得出四个点偏差太大，可以考虑去除这四个点
 
subplot(211);                                                           % 画残差图
rcoplot(r,rint);                                                        % 画残差图
 
第四步：得出预测模型，用预测模型去和真实值去对比
 
 
NewTotalEnergy=zeros(Row,1);                                            % 预测数据矩阵
for row = 1:Row                                                         % 预测赋值
    NewTotalEnergy(row)=b(1);
    for PaRow=2:(count+1)
        NewTotalEnergy(row)=NewTotalEnergy(row)+b(PaRow)*ConVariable(row,PaRow-1);
    end
end
Loss = 0;                                   % 根据最大斜率和最小效率建立指定长度的损失数据       
for i=1:Row
    Loss = Loss+(TotalEnergy(i)-NewTotalEnergy(i))^2/(2*Row);
end
subplot(212);                                                         % 画预测和真实图
number=[1:1:Row];
plot(number,TotalEnergy','r',number,NewTotalEnergy','b');
xlabel('数量序列');ylabel('总能耗');title('总能耗与相关参数散点图');legend('真实值','预测值')grid on;

 
第五步：得出模型函数关系
 $$总能耗=33374+-53041\times 体型系数+32\times 面积-101\times 人口密度+431\times 内扰电耗$$
 
参考文献
 
数学建模与数学试验
多元线性回归MATLAB实现
一元线性回归模型及其假设条件
[MATLAB]逐步回归详解(stepwise使用指南)
基于Matlab的数据多元回归分析的研究

1、案例
 
%%1、bint表示回归系数区间估计可参考http://www.360doc.com/content/11/0801/20/2537127_137246007.shtml
%2、r表示残差
%3、rint代表置信区间
%4、stas表示用于检验回归模型的统计量，有三个数值 r^2 F 与F对应的概率P 例如p<0.05 残差95%
%   r^2越接近于1，回归方程越显著  
%alpha表示显著水平

%%
x=[143 144 145 147 148 150 153 154 155 156 157 158 159 160 161 162]';
X=[ones(16,1),x];
Y=[87 85 88 91 92 90 93 95 98 98 97 95 97 99 100 102]';
[b,bint,r,rint,stats]=regress(Y,X)
t=1:16;
%%
figure(1);
y_fitting=X(t,:)*b;
plot(t,y_fitting,'r-',  t,Y(t,:),'b-', t,abs(y_fitting-Y(t,:)),'k-');
legend('红--拟合值','蓝--实际值','黑--误差值');
text(3,50,strcat('相关系数R=',num2str(stats(1,1 ))));
text(7,50,strcat('F=',num2str(stats(1,2))));
text(9,50,strcat('P=',num2str(stats(1,3 ))));
nhfcs1=strcat('拟合方程式',num2str(b(1,1)),'+',num2str(b(2,1)),'*X1');
text(11,50,nhfcs1);
%
%功能 在当前轴中创建text对象。函数text是创建text图形句柄的低级函数。可用该函数在图形中指定的位置上显示字符串。

%用法 text(x,y,'string')在图形中指定的位置(x,y)上显示字符串string

%text(x,y,z,'string') 在三维图形空间中的指定位置(x,y,z)上显示字符串string
%
title('线性回归曲线拟合结果');
xlabel('样本点');
ylabel('分数');

%%
figure(2);
ul=rint(:,1);
I1=rint(:,2);
plot(t,I1,'b-', t,r,'R*',  t, ul,'g-');
legend('蓝色--残差95%置信区间上限','红--残差值','绿--残差95%置信区间下限');
xlabel('样本值');
ylabel('残差值');
figure(3)
rcoplot(r,rint);   %残差分析，作残差图
 

 
 
 
 
 
2、残差图(Residual Plots)
 
我们可以用残差图来估计观察或预测到的误差error(残差residuals)与随机误差(stochastic error)是否一致。用一个丢骰子的例子最好理解了。当你丢出去一个六面的骰子时，你不应该能够预测得到哪面点数向上。然而，你却可以评估在一系列投掷后，正面向上的数字是否遵循一个随机模式，你自己心中就会想象出一个随机散布的残差图。如果，有人背着你对骰子做了点手脚，让六点更频繁的出现向上，这时你心中的残差图看上去就似乎有规律可循，从而不得不修改心中的模型，让你狐疑骰子一定有问题。
 相同的原则也适用于回归模型。你不应该能够预测任何给定的观察或预测结果的错误(或者说差别)。你需要确定残差是否与随机误差相互呈现一致性，就像丢骰子一样，残差若整体呈现“很古怪”的模式，你就需要回头修改你的回归模型了。上面“古怪”究竟怎么看呢？看下文。
 
 
%clc
%clear
%%
%目标函数：y=Ax1^2+Bx2^2+Cx1+Dx2+Ex1*x2+F  （这是一个二次函数，两个变量，大写的字母是常数）
%导入数据  
y=[7613.51  7850.91  8381.86  9142.81 10813.6 8631.43 8124.94 9429.79 10230.81 10163.61 9737.56 8561.06 7781.82 7110.97]';  
x1=[7666 7704 8148 8571 8679 7704 6471 5870 5289 3815 3335 2927 2758 2591]';  
x2=[16.22 16.85 17.93 17.28 17.23 17 19 18.22 16.3 13.37 11.62 10.36 9.83 9.25]';  
X=[ones(size(y)) x1.^2 x2.^2 x1 x2 x1.*x2];  
%开始分析  
[b,bint,r,rint,stats] = regress(y,X)
scatter3(x1,x2,y,'filled') %scatter可用于画散点图a
hold on
%%
%拟合，三维视图显示  
hold on  %不要清除计算数据，在刚刚那副散点图上接着画  
x1fit = min(x1):100:max(x1);   %设置x1的数据间隔  
x2fit = min(x2):1:max(x2);     %设置x2的数据间隔  
[X1FIT,X2FIT] = meshgrid(x1fit,x2fit);  %生成一个二维网格平面，也可以说生成X1FIT,X2FIT的坐标  
YFIT=b(1)+b(2)*X1FIT.^2+b(3)*X2FIT.^2+b(4)*X1FIT+b(5)*X2FIT+b(6)*X1FIT.*X2FIT;    %代入已经求得的参数，拟合函数式  
mesh(X1FIT,X2FIT,YFIT)    %X1FIT，X2FIT是网格坐标矩阵，YFIT是网格点上的高度矩阵  
view(10,10)  %改变角度观看已存在的三维图，第一个10表示方位角，第二个表示俯视角。  
             %方位角相当于球坐标中的经度，俯视角相当于球坐标中的纬度  
xlabel('x1') %设置X轴的名称  
ylabel('x2') %设置y轴的名称  
zlabel('y')  %设置z轴的名称
hold on



%%
figure(2)
rcoplot(r,rint);   %残差分析，作残差图
 

 
 
 
 
https://zhuanlan.zhihu.com/p/20700731
 
3、回归值与残差的残差图编辑
 
为检验建立的多元线性回归模型是否合适，可以通过回归值与残差的散点图来检验。其方法是画出回归值与普通残差的散点图，或者画出回归值与标准残差的散点图，其图形可能会出现下面三种情况(如图1所示)：
 
图1（a）
 
图1（b）
 
对于图1(a)的情况，不论回归值的大小，而残差(或)具有相同的分布，并满足模型的各假设条件；对于图1(b)的情况，表示回归值的大小与残差的波动大小有关系，即等方差性的假设有问题；对于图1©，表示线性模型不合适的样本，可能有异常值存在。
 
对于图1(a)，如果大部分点都落在中间(b)部分，而只有少数几个点落在外边，则这些点对应的样本，可能有异常值存在。
 
图1（c）