实验11_统计回归模型(4学时)

1实验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.3823.754.006.750.258.5133.704.307.250.609.5243.703.705.5007.5053.603.857.000.259.3363.603.806.500.208.2873.603.756.750.158.7583.803.855.250.057.8793.803.655.25-0.157.10103.854.006.000.158.00113.904.106.500.207.89123.904.006.250.108.15《数学建模实验》王平2133.704.107.000.409.10143.754.206.900.458.86153.754.106.800.358.90163.804.106.800.308.87173.704.207.100.509.26183.804.307.000.509.00193.704.106.800.408.75203.803.756.50-0.057.95213.803.756.25-0.057.65223.753.656.00-0.107.27233.703.906.500.208.00243.553.657.000.108.50253.604.106.800.508.75263.654.256.800.609.21273.703.656.50-0.058.27283.753.755.7507.67293.803.855.800.057.93303.704.256.800.559.261.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。3☆(1)运行程序并给出结果（比较[327]图1）：(2)确定y对x1的拟合，绘制散点图与拟合曲线组合图形从y对x1的散点图可以发现，可用线性模型（直线）011yx来拟合（其中ε是随机误差）。程序如下：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;4[提示：regress多元线性回归函数调用格式][b,bint,r,rint,stats]=regress(y,x,alpha)例，多元回归模型为：20112232yxxx输入：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剩余方差5☆(2)运行程序并给出结果（比较[327]图1）：(3)绘制y对x2的散点图程序如下：clc;formatshortg;M=dlmread('p325.txt');%读取ASCII码文件x2=M(:,4);y=M(:,6);plot(x2,y,'bo');6☆(3)运行程序并给出结果（比较[327]图2）：(4)确定y对x2的的拟合，绘制散点图与拟合曲线组合图形从y对x2的散点图可以发现，可用二次函数模型201222yxx来拟合。程序如下：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;7☆(4)运行程序并给出结果（比较[327]图2）：(5)y对x1,x2的回归模型及其求解，销售量预测综上得回归模型20112232yxxx变量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);8x1=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的回归模型的改进和求解，销售量预测改进的模型20112232412yxxxxx参考题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位小数9fprintf('\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)（验证）完全二次多项式模型22011223124152yxxxxxx运行以下程序（参考教材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估计值和预测区间。10点击所得交互画面左下方的输出按钮“Export”，所得画面（导出到工作空间）第1个复选框是“将拟合参数存到一个名为beta的MATLAB变量中”，点击OK。在命令窗口提示符键入变量名beta将得到参数(β0,β1,β2,β3,β4,β5)'的值。☆(2)运行程序并给出结果：参数(β0,β1,β2,β3,β4,β5)'的值（比较[331]）：112.软件开发人员的薪金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）01138761110211608103031870111304112831020511767103062087221207117722020810535201091219520310123133021114975311122137131213198003131411417401152026341316132314031712884402181324550219136775031220159655112112366601222135261323138396022422884612251697871126148038022717404811282218481329135488013014467100131159421002322317410133323780101234254101112351486111013616882120237241701213381599013013926330131240179491402412568515134227837161243188381602441748316014519207170246193462001假设满足多元线性回归模型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','参数','估计值','置信区间'