实验11统计回归模型4学时.docx
《实验11统计回归模型4学时.docx》由会员分享，可在线阅读，更多相关《实验11统计回归模型4学时.docx（68页珍藏版）》请在冰豆网上搜索。
实验11统计回归模型4学时
实验11统计回归模型（4学时）
（第10章统计回归模型）
1.牙膏的销售量p325~332
下面给出一组数据，其中：
第1列销售周期；
第2列某公司牙膏销售价格（元）x4；
第3列其它厂家平均价格（元）x3；
第4列广告费用（百万元）x2；
第5列价格差（元）x1（x3-x4）；
第6列销售量（百万支）y。
存放在一个名为p325.txt的文件中。
13.853.805.50-0.057.38
23.754.006.750.258.51
33.704.307.250.609.52
43.703.705.5007.50
53.603.857.000.259.33
63.603.806.500.208.28
73.603.756.750.158.75
83.803.855.250.057.87
93.803.655.25-0.157.10
103.854.006.000.158.00
113.904.106.500.207.89
123.904.006.250.108.15
133.704.107.000.409.10
143.754.206.900.458.86
153.754.106.800.358.90
163.804.106.800.308.87
173.704.207.100.509.26
183.804.307.000.509.00
193.704.106.800.408.75
203.803.756.50-0.057.95
213.803.756.25-0.057.65
223.753.656.00-0.107.27
233.703.906.500.208.00
243.553.657.000.108.50
253.604.106.800.508.75
263.654.256.800.609.21
273.703.656.50-0.058.27
283.753.755.7507.67
293.803.855.800.057.93
303.704.256.800.559.26
1.1（验证）基本模型p325~329
先保存上面的p325.txt文件。
（1）绘制y对x1的散点图
程序如下：
M=dlmread（'p325.txt'）;%读取ASCII码文件
x1=M（:
5）;y=M（:
6）;
plot（x1,y,'bo'）;
[提示：
dlmread将以ASCII码分隔的数值数据文件读入到矩阵中]
dlmread：
读取ASCII码文件的MATLAB函数
M=dlmread（'fun.txt'）;
fun.m是一个数据文件，存放一个数据矩阵，将文件内容写入M。
☆
（1）运行程序并给出结果（比较[327]图1）：
（2）确定y对x1的拟合，绘制散点图与拟合曲线组合图形
从y对x1的散点图可以发现，可用线性模型（直线）
来拟合（其中ε是随机误差）。
程序如下：
clc;formatshortg;
M=dlmread（'p325.txt'）;%读取ASCII码文件
x1=M（:
5）;y=M（:
6）;
plot（x1,y,'bo'）;
b=regress（y,[ones（size（x1））,x1]）;%b=[β0β1]'，列向量
x1=sort（x1）;%按升序排序，用于画图
y=[ones（size（x1））,x1]*b;%使用矩阵乘法
holdon;
plot（x1,y,'-r'）;
holdoff;
[提示：
regress多元线性回归函数调用格式]
[b,bint,r,rint,stats]=regress（y,x,alpha）
例，多元回归模型为：
输入：
y为n（=30）维列向量数据。
x为对应于回归系数（β0,β1,β2,β3）'的数据矩阵[1x1x2x22]（30×4矩阵，第1列全1）。
alpha为置信水平（缺省时为0.05）。
输出：
b为β=（β0,β1,β2,β3）'估计值，4维列向量。
bint为b的置信区间，4×2矩阵。
r为残差n（=30）维列向量y-xβ。
rint为r的置信区间，30×2矩阵。
stats为回归模型的检验统计量，含4个值：
R2回归方程的决定系数（R是相关系数）
F统计值
P与F统计量对应的概率值
s2剩余方差
☆
（2）运行程序并给出结果（比较[327]图1）：
（3）绘制y对x2的散点图
程序如下：
clc;formatshortg;
M=dlmread（'p325.txt'）;%读取ASCII码文件
x2=M（:
4）;y=M（:
6）;
plot（x2,y,'bo'）;
☆（3）运行程序并给出结果（比较[327]图2）：
（4）确定y对x2的的拟合，绘制散点图与拟合曲线组合图形
从y对x2的散点图可以发现，可用二次函数模型
来拟合。
程序如下：
clc;formatshortg;
M=dlmread（'p325.txt'）;%读取ASCII码文件
x2=M（:
4）;y=M（:
6）;
plot（x2,y,'bo'）;
b=regress（y,[ones（size（x2））,x2,x2.^2]）;%b=[β0β1β2]'，列向量
x2=sort（x2）;
y=[ones（size（x2））,x2,x2.^2]*b;%使用矩阵乘法
holdon;
plot（x2,y,'-r'）;
holdoff;
☆（4）运行程序并给出结果（比较[327]图2）：
（5）y对x1,x2的回归模型及其求解，销售量预测
综上得回归模型
变量x1,x2为回归变量，参数β0,β1,β2,β3为回归系数。
程序如下：
clc;formatcompact;formatshortg;
M=dlmread（'p325.txt'）;%读取ASCII码文件
x1=M（:
5）;x2=M（:
4）;y=M（:
6）;
[b,bint,r,rint,stats]=regress（y,[ones（size（x1））,x1,x2,x2.^2],0.05）;
fprintf（'%2s%5s%11s\n','参数','估计值','置信区间'）;%1个汉字算1个字符
fori=1:
length（b）
fprintf（'β%1d%9.4f[%7.4f,%7.4f]\n',i-1,[b（i,:
）,bint（i,:
）]）;
end%%d将i当整数输出，%7.4f按实数格式输出，区域宽7个字符，4位小数
fprintf（'\nR2=%.4fF=%.4fp<%.4es2=%.4f\n',stats）;
x1=0.2;x2=6.5;y=[1x1,x2,x2^2]*b;%使用矩阵乘法
fprintf（'\n销售量预测：
x1=%.1f,x2=%.1f,y=%.4f\n',x1,x2,y）;
[提示：
fprintf输出到命令窗口或写数据到文本文件]
见参考资料：
MATLAB函数和命令的用法。
☆（5）运行程序并给出结果（比较[328]表2,[329]的预测结果）：
1.2（验证，编程）模型改进p329~332
仍使用题1的数据。
（1）（编程）y对x1,x2的回归模型的改进和求解，销售量预测
改进的模型
参考题1（5）的程序，编写一个类似的程序，运行结果与教材p329~330的表3及相关结果相比较。
★
（1）给出程序和运行结果（比较[329]表3）：
clc;formatcompact;formatshortg;
M=dlmread（'p325.txt'）;%读取ASCII码文件
x1=M（:
5）;x2=M（:
4）;y=M（:
6）;
[b,bint,r,rint,stats]=regress（y,[ones（size（x1））,x1,x2,x2.^2,x1.*x2],0.05）;
fprintf（'%2s%5s%11s\n','参数','估计值','置信区间'）;%1个汉字算1个字符
fori=1:
length（b）
fprintf（'β%1d%9.4f[%7.4f,%7.4f]\n',i-1,[b（i,:
）,bint（i,:
）]）;
end%%d将i当整数输出，%7.4f按实数格式输出，区域宽7个字符，4位小数
fprintf（'\nR2=%.4fF=%.4fp<%.4es2=%.4f\n',stats）;
x1=0.2;x2=6.5;y=[1x1,x2,x2^2,x1.*x2]*b;%使用矩阵乘法
fprintf（'\n销售量预测：
x1=%.1f,x2=%.1f,y=%.4f\n',x1,x2,y）;
（2）（验证）完全二次多项式模型
运行以下程序（参考教材p331~332）：
clear;clc;formatcompact;formatshortg;
M=dlmread（'p325.txt'）;%读取ASCII码文件
x1=M（:
5）;x2=M（:
4）;y=M（:
6）;
rstool（[x1,x2],y,'quadratic'）
得以下的交互画面。
画面中的两个座标系给出y的估计值和预测区间。
用鼠标移动交互式画面中的十字线，或在图下方的窗口内输入，可改变x1和x2的数值。
改变x1=0.2，x2=6.5，观察窗口左边的y估计值和预测区间。
点击所得交互画面左下方的输出按钮“Export”，所得画面（导出到工作空间）第1个复选框是“将拟合参数存到一个名为beta的MATLAB变量中”，点击OK。
在命令窗口提示符键入变量名beta将得到参数（β0,β1,β2,β3,β4,β5）'的值。
☆
（2）运行程序并给出结果：
参数（β0,β1,β2,β3,β4,β5）'的值（比较[331]）：
2.软件开发人员的薪金p332~338
在下面给出的数据中：
（存入文件p333.txt）
第1列编号
第2列薪金y
第3列资历x1（从事专业工作的年数）
第4列管理x2（1表示管理人员，0表示非管理人员）
第5列教育x3，x4（1表示中学程度x3x4=10，2为大学x3x4=01，3为更高程度x3x4=00）
0113876111
0211608103
0318701113
0411283102
0511767103
0620872212
0711772202
0810535201
0912195203
1012313302
1114975311
1221371312
1319800313
1411417401
1520263413
1613231403
1712884402
1813245502
1913677503
2015965511
2112366601
2221352613
2313839602
2422884612
2516978711
2614803802
2717404811
2822184813
2913548801
30144671001
31159421002
32231741013
33237801012
34254101112
35148611101
36168821202
37241701213
38159901301
39263301312
40179491402
41256851513
42278371612
43188381602
44174831601
45192071702
46193462001
假设满足多元线性回归模型
y=α0+α1x1+α2x2+α3x3+α4x4+ε
2.1（验证）基本模型p332~335
求回归系数及其置信区间（置信水平α=0.05）、检验统计量R2、F、p、s2，有关散点图的MATLAB程序如下：
%10.2软件开发人员的薪金——基本模型
%模型：
y=α0+α1*x1+α2*x2+α3*x3+α4*x4+ε
clear;clc;formatcompact;formatshortg;
M=dlmread（'p333.txt'）;%读取ASCII码文件
y=M（:
2）;x1=M（:
3）;x2=M（:
4）;
x3=M（:
5）==1;x4=M（:
5）==2;%教育程度
[a,aint,r,rint,stats]=regress（y,[ones（size（M,1）,1）x1x2x3x4]）;
fprintf（'%2s%4s%9s\n','参数','估计值','置信区间'）;%1个汉字算1个字符
fori=1:
length（a）
fprintf（'α%1d%7.0f[%5.0f,%5.0f]\n',i-1,[a（i,:
）,aint（i,:
）]）;
end
fprintf（'\nR2=%.3fF=%.0fp<%.4es2=%.3e\n',stats）;
subplot（121）;
plot（x1,r,'+'）;%模型的残差ε与资历x1的关系
subplot（122）;
x234=（1+M（:
4））.*（M（:
5）==1）+（3+M（:
4））.*（M（:
5）==2）...
+（5+M（:
4））.*（M（:
5）==3）;%见p336表3
plot（x234,r,'+'）;%模型的残差ε与管理-教育x2-x3,x4组合x234的关系
☆给出程序的运行结果（数据和图形）（比较[334]，[335]）：
数据结果（比较[334]表2）：
图形结果（比较[335]图1、图2）：
2.2（编程）更好的模型p335~336
在题2.1的模型中增加x2与x3，x4的交互项后，多元线性回归模型为
y=α0+α1x1+α2x2+α3x3+α4x4+α5x2x3+α6x2x4+ε
要求：
同题2.1，区别在于模型不同，所以要根据新模型修改题2.1的程序，仍后运行。
★给出程序和运行结果（程序、数据和图形）（比较[336]表4、图3、图4）：
程序：
%10.2软件开发人员的薪金——基本模型
%模型：
y=α0+α1*x1+α2*x2+α3*x3+α4*x4+ε
clear;clc;formatcompact;formatshortg;
M=dlmread（'p333.txt'）;%读取ASCII码文件
y=M（:
2）;x1=M（:
3）;x2=M（:
4）;
x3=M（:
5）==1;x4=M（:
5）==2;%教育程度
[a,aint,r,rint,stats]=regress（y,[ones（size（M,1）,1）x1x2x3x4x2.*x3x2.*x4]）;
fprintf（'%2s%4s%9s\n','参数','估计值','置信区间'）;%1个汉字算1个字符
fori=1:
length（a）
fprintf（'α%1d%7.0f[%5.0f,%5.0f]\n',i-1,[a（i,:
）,aint（i,:
）]）;
end
fprintf（'\nR2=%.3fF=%.0fp<%.4es2=%.3e\n',stats）;
subplot（121）;
plot（x1,r,'+'）;%模型的残差ε与资历x1的关系
subplot（122）;
x234=（1+M（:
4））.*（M（:
5）==1）+（3+M（:
4））.*（M（:
5）==2）...
+（5+M（:
4））.*（M（:
5）==3）;%见p336表3
plot（x234,r,'+'）;%模型的残差ε与管理-教育x2-x3,x4组合x234的关系
数据结果（比较p336表4）：
图形结果（比较p336图3、图4）：
2.3（编程，验证）模型应用p337~338
在题2.1的模型中增加x2与x3，x4的交互项后，多元线性回归模型为
y=α0+α1x1+α2x2+α3x3+α4x4+α5x2x3+α6x2x4+ε
（1）（编程）同题2.2，但要去除题2.1中给出数据中编号为33的数据（异常点）
★
（1）给出程序和运行结果（数据和图形）（比较[337]表5、图5、图6）：
程序：
%10.2软件开发人员的薪金——基本模型
%模型：
y=α0+α1*x1+α2*x2+α3*x3+α4*x4+ε
clear;clc;formatcompact;formatshortg;
M=dlmread（'p333.txt'）;%读取ASCII码文件
M（33,:
）=[];
y=M（:
2）;x1=M（:
3）;x2=M（:
4）;
x3=M（:
5）==1;x4=M（:
5）==2;%教育程度
[a,aint,r,rint,stats]=regress（y,[ones（size（M,1）,1）x1x2x3x4x2.*x3x2.*x4]）;
fprintf（'%2s%4s%9s\n','参数','估计值','置信区间'）;%1个汉字算1个字符
fori=1:
length（a）
fprintf（'α%1d%7.0f[%5.0f,%5.0f]\n',i-1,[a（i,:
）,aint（i,:
）]）;
end
fprintf（'\nR2=%.3fF=%.0fp<%.4es2=%.3e\n',stats）;
subplot（121）;
plot（x1,r,'+'）;%模型的残差ε与资历x1的关系
subplot（122）;
x234=（1+M（:
4））.*（M（:
5）==1）+（3+M（:
4））.*（M（:
5）==2）...
+（5+M（:
4））.*（M（:
5）==3）;%见p336表3
plot（x234,r,'+'）;%模型的残差ε与管理-教育x2-x3,x4组合x234的关系
数据结果（比较p337表5）：
图形结果（比较p337图5、图6）：
（2）（验证）6种管理—教育组合人员的“基础”薪金（资历为0）。
组合
管理教育“基础”薪金
资历
1
01？
0
2
11？
0
3
02？
0
4
12？
0
5
03？
0
6
13？
0
利用
（1）所求得的模型求解，程序如下（与教材p337表6比较）：
clear;clc;formatcompact;formatshortg;
M=dlmread（'p333.txt'）;%读取ASCII码文件
M（33,:
）=[];%修改
y=M（:
2）;x1=M（:
3）;x2=M（:
4）;
x3=M（:
5）==1;x4=M（:
5）==2;%教育程度
[a,aint,r,rint,stats]=regress（y,[ones（size（M,1）,1）x1x2x3x4x2.*x3x2.*x4]）;
x1=zeros（6,1）;%资历为0
x2=[010101]';%管理
x3=[110000]';%教育
x4=[001100]';
y=[ones（6,1）,x1,x2,x3,x4,x2.*x3,x2.*x4]*a;
fprintf（'%2s%4s%4s%8s\n','组合','管理','教育','“基础”薪金'）;%1个汉字算1个字符
fork=1:
6
fprintf（'%3d%6d%6d%11.0f\n',k,x2（k）,fix（（k+1）/2）,y（k））;
end
☆
（2）运行程序并给出结果（比较[337]表6）：
3.酶促反应p338~345
嘌呤霉素实验中的反应速度与底物浓度数据
底物浓度
（ppm）
0.02
0.06
0.11
0.22
0.56
1.10
反应
速度
处理
76
47
97
107
123
139
159
152
191
201
207
200
未处理
67
51
84
86
98
115
131
124
144
158
160
/
3.1（编程）线性化模型p338~340
（1）根据给出的数据，试编写一个程序绘制出如下图所示的散点图（比较：
p339图1、图2）
★
（1）给出程序和运行结果（比较[339]图1、图2）：
程序：
clear;clc;formatcompact;formatshortg;
M=dlmread（'p338.txt'）;%读取ASCII码文件
y=M（2,1:
2:
11）;y2=M（2,2:
2:
12）;x1=M（1,1:
6）;x2=M（1,1:
5）;
y4=M（3,1:
2:
11）;y3=M（3,2:
2:
10）;
subplot（1,2,1）
plot（x1,y,'bo'）;
holdon;
plot（x1,y2,'bo'）;
title（'y对x经处理的散点图'）;
subplot（1,2,2）
plot（x1,y4,'+'）;
holdon;
plot（x2,y3,'+'）;
title（'y对x未经处理的散点图'）;
图形（比较p339图1、图2）：
（2）线性化模型
Michaelis-Menten模型
通过代换
或
可化为线性模型
试编写一个程序，对经过处理的实验数据，求出：
（a）参数的估计；
（b）1/y与1/x的散点图和回归直线；
（c）用线性化得到的原始数据拟合图。
程序运行的结果如下：
p340表2
p339图3、p340图4
★
（2）给出程序和运行结果（比较[339]，[340]）：
程序：
数据结果（比较p340表2）：
图形（比较p339图3、p340图4）：
3.2（验证）非线性模型及求解p340~341
Michaelis-Menten模型
[提
升级会员 