系统工程导论章四作业一

《系统工程导论》第四章作业

利用Matlab，编程实现一元线性回归。

要求：

1）实现函数function linear_regression1(data,alpha);

2）输入为N*2的矩阵data，第一列为Y，第二列为X;显著性水平alpha；

3）打印出回归直线方程（有必要的话，也可输出重要的中间数据）；

4）用F检验法进行统计检验，显著性水平为输入alpha，提示Matlab中得到F分布对于给定显著性水平和自由度的分位数函数是finv，请大家自行用help工具确定其用法；输出检验结果，如果输入数据满足线性关系，那么5）和6），否则结束；

5）打印出置信区间，提示Matlab中得到标准正态分布相应分位数的函数是norminv，请大家自行用help工具确定其用法；

6）在一个figure中，画出:a所有数据点，b回归直线，c置信区间相应的两条边界直线。

7）将“第4章作业(1)_数据.ppt”中的数据用你

编写的程序处理，将所有结果贴到作业报告中，

显著性水平取0.05。

解：

算例结果如下图：

打印出回归直线

用F检验法进行统计检验，并输出检验结果

打印出置信区间

在一个figure中，画出:a所有数据点，b回归直线，c置信区间相应的两条边界直线：

注：图中，红点是数据点，红线是回归直线，蓝线是置信区间相应的两条边界直线。

代码：

function linear_regressionl(data,alpha);

data=[4.0 3.44 3.6 1.0 2.04 4.74 0.6 1.7 2.92 4.8 3.28 4.16 3.35 2.2;

0.009 0.013 0.006 0.025 0.022 0.007 0.036 0.014 0.016 0.014 0.016 0.012

0.020 0.018];//输入N*2的矩阵data，第一列为Y第二列为X

alpha=0.05; //置信水平

X=[ 1 1 1 1 1 1 1 1 1 1 1

1 1 1;

0.009 0.013 0.006 0.025 0.022 0.007 0.036 0.014 0.016 0.014 0.016 0.012 0.020 0.018];

Y=[4.0 3.44 3.6 1.0 2.04 4.74 0.6 1.7 2.92 4.8 3.28 4.16 3.35 2.2];

A=X*X';

B=X*Y';

M=inv(A)*B; //M是回归模型式中的参数

disp('y=5.18-134.6066*x'); //输出回归直线

Fa=finv(0.95,1,12) //求出临界值

ESS=0;RSS=0;yi=0;xi=0;

for i=1:14

yi=yi+Y(i);

xi=xi+X(2,i);

end

y=yi/14;

x=xi/14; //分别求出x与y的平均值

for i=1:14

ESS=ESS+((5.18-134.6066*data(2,i))-y)*((5.18-134.6066*data(2,i))-y);

RSS=RSS+(Y(i)-(5.18-134.6066*data(2,i)))*(Y(i)-(5.18-134.6066*data(2,i))); end

F=12*ESS/RSS

disp('F>Fa 所以是线性的');

Z=norminv(0.025) //标准正态分布相应分位数

S=sqrt(RSS/12) //剩余均方差

disp('置信区间：（y0-1.59，y0+1.59）');

x=0.005:0.0001:0.04; //取x的范围和步长

y=5.18-134.6066*x;

m=[0.009 0.013 0.006 0.025 0.022 0.007 0.036 0.014 0.016 0.014 0.016 0.012 0.020 0.018];

n=[4.0 3.44 3.6 1.0 2.04 4.74 0.6 1.7 2.92 4.8 3.28 4.16 3.35 2.2];

//为了方便输出数据点方便，用m与n分别表示x，y。因为x、y存储在向量里面不好用

plot(m,n,'r.'); //输出数据点

hold on;

plot(x,y,'r'); //输出回归直线

hold on;

plot(x,y+Z*S,'b');

hold on;

plot(x,y-Z*S,'b'); //输出置信区间相应的两条边界直线。

for i=1:14

str=['(',num2str(data(2,i)),',',num2str(data(1,i)),')'];

text(data(2,i),data(1,i),str); //输出坐标

end

title('y=5.18-134.6066*x'); //写标题

xlabel ('x');

ylabel ('y'); //标注横坐标、纵坐标

end