数学建模举例

10.1牙膏的销售量
某大型牙膏制造企业为了更好地拓展产品市场，有效地管理库存，公司董事会要求销售部门根据市场调查，找出公司生产的牙膏销售量与销售价格、广告投入等之间的关系，从而预测出在不同价格和广告费用下的销售量。

为此，销售部的研究人员收集了过去30个销售周期（每个销售周期为4周）公司生产的牙膏的销售量、销售价格、投入的广告费用，以及同期其它厂家生产的同类牙膏的市场平均销售价格，见表1-1（其中价格差指其它厂家平均价格与公司销售价格之差）。

试根据这些数据建立一个数学模型，分析牙膏销售量与其它因素的关系，为制订价格策略和广告投入策略提供数量依据
表1-1牙膏销售量与销售价格、广告费用等数据
一、问题重述
根据过去30个销售周期（每个销售周期为4周）公司生产的牙膏的销售量、销售价格、投入的广告费用，以及同期其它厂家生产的同类牙膏的市场平均销售价格，见表1-1。

根据这些数据建立一个数学模型，分析牙膏销售量与其它因素的关系，为制订价格策略和广告投入策略提供数量依据
二、问题分析
由于牙膏是生活必需品，对大多属顾客来说，在购买同类产品的牙膏是更多地会在意不同品牌之间的价格差异，而不是它们的价格本身。

因此，在研究各个因素对销量的影响时，用价格差代替公司销售价格和其他厂家平均价格更为合适。

三、模型假设
1.画出牙膏销售量与价格差，公司投入的广告费用的散点图
2.由散点图确定两个函数模型，再由这两个函数模型解出回归模型
3.对模型进行改进，添加新的条件确定更好的回归模型系数，得到新的回归模型
4.对模型进一步改进，确定最终的模型
四、符号约定
牙膏销售量为y，其他厂家平均价格和公司销售价格之差（价格差）为x1，公司投入的广告费用为x2，其他厂家平均价格和公司销售价格分别为x3和x4,x1=x3-x4。

基于上面的分析，我们仅利用1x和2x来建立y的预测模型。

五、模型的建立和求解
1.基本模型
利用表1-1的数据用matlab作出y与x1的散点图（图1-1），y与x2的散点图（图1-2）
代码如下：
x1=[-0.05 0.25 0.6 0 0.25 0.2 0.15 0.05 -0.15 0.15 0.2 0.1 0.4 0.45 0.35 0.3 0.5 0.5 0.4 -0.05 -0.05 -0.1 0.2 0.1 0.5 0.6 -0.05 0 0.05 0.55];
x2=[5.5 6.75 7.25 5.5 7 6.5 6.75 5.25 5.25 6 6.5 6.25 7 6.9 6.8 6.8 7.1 7 6.8 6.5 6.25 6 6.5 7 6.8 6.8
6.5 5.75 5.8 6.8];
y=[7.38 8.51 9.52 7.5 9.33 8.28 8.75 7.87 7.1 8 7.89 8.15 9.1 8.86 8.9 8.87 9.26 9 8.75 7.95 7.65
7.27 8 8.5 8.75 9.21 8.27 7.67 7.93 9.26];
A1=polyfit(x1,y,1);
yy1=polyval(A1,x1);
A2=polyfit(x2,y,2);
x5=5:0.05:7.25;
yy2=polyval(A2,x5);
subplot(1,2,1);plot(x1,y,'o',x1,yy1);
title('图1 y对x1的散点图');
subplot(1,2,2);plot(x2,y,'o',x5,yy2);
title('图2 y对x2的散点图');
图（1-1）与图（1-2）
从图1可以发现，随着1x 的增加，y 的值有比较明显的线性增长趋势，图中的直线是用线性模型：
011y x ββε=++（1）
拟合的（其中ε是随机误差）。

而在图2中，当x2增大时，y 有向上弯曲增加的趋势，图中的曲线是用二次函数模型:
021222y x x βββε=+++（2）
综合上面的分析，结合模型（1）和（2）建立如下的回归模型:
02112232y x x x ββββε=++++（3）
（3）式右端的x1和x2称为回归变量（自变量），02
112232x x x ββββ+++是给定价格差x1，广告费用x2时，牙膏销售量y 的平均值，其中的参数0123,,,ββββ称为回归系数，由表1-1的数据估计，影响y 的其他因素作用都包含在随机误差ε中。

如果模型选择合适，ε应该大致服从值为0的正态分布。

2. 模型求解
在刚刚运行的代码后面，继续使用regress 工具求解，代码为：
x6=[ones(30,1) x1' x2' (x2.^2)']; [b,bint,r,rint,stats]=regress(y',x6,0.05) 运行结果如图（1-3）
得到模型（3）的回归系数估计值及其置信区间（置信水平α=0.05）、检验统计量2
R ,F ，p ，2s 得结果见表1-2
图（1-3）
3. 结果分析
表1-2显示，2
R =0.9054指因变量y （销售量）的90.54%可由模型决定，F 值远远超过F 检验的临界值，P 远小于α，因而模型（3）整体来看是可用的
表1-2的回归系数中
2β的置信区间包含零点，
表示回归变量2
x （对因变量y 的影响）是不太显著的，
但由于2
x 是显著的，我们仍将2
x 留在模型中
4. 销售量预测
将回归系数的估计值带入模型（3），即可预测公司未来某个销售周期牙膏的销售量y ，预测值记为y ∧
，得到模型（3）的预测方程：
02
123122y x x x ββββ∧
∧
∧
∧
=+++（4）
只需要知道该销售周期的价格差x1和投入的广告费用x2,就可以计算预测值y ∧。

5. 模型改进
模型（3）中回归变量x1和x2对因变量y 的影响是相互独立的，即牙膏销售量y 的均值与广告费用
x2的二次关系由回归系数2β和3β确定，而不依赖于价格差x1，同样的，y 的均值与x1的线性关系由回归
系数
1β确定，而不依赖于x2。

根据直觉和经验可以猜想，x1和x2之间的交互作用会对y 有影响，不妨简
单地用x1，x2的乘积代表它们的相互作用，于是将模型（3）增加一项，得到:
021********y x x x x x βββββε=+++++（5）
在这个模型中，y 的均值与2x 的二次关系为
2241232x x x βββ++（），由系数确 定，并依赖于价格差x1。

在上述运行程序后继续输入代码：
x7=[ones(30,1) x1' x2' (x2.^2)' (x1.*x2)']; [b,bint,r,rint,stats]=regress(y',x7,0.05); b,bint,stats 结果见图（1-4）
图（1-4）
R有所提高，说明模型（5）比模型（3）有所进步。

并且，所有参数的置表3与表2的结果相比，2
的置信区间不包含零点，所以有理由相信模型（5）信区间，特别是X1，X2的交互作用项X1X2的系数
4
比模型（3）更符合实际。

在保持广告费用x2=6.5百万元不变的条件下，分别对模型（3）和（5）中牙膏销售量的均值与价格差x1的关系作图，见图1-5和图1-6，代码为：
yy3=17.3244+1.307*x1+(-3.6956)*6.5+0.3486*6.5*6.5;
plot(x1,yy3);
grid on
figure(2)
yy4=29.1133+11.1342*x1+(-7.608*6.5)+0.6712*6.5*6.5+(-1.4777)*6.5*x1;
plot(x1,yy4);
grid on
图1-5
图1-6
在保持价格差x1=0.2元不变的条件下，分别对模型（3）和（5）中牙膏销售量的均值y
与广告费用x2的关系作图，见图1-7和图1-8，代码如下：
figure(3)
yy5=17.3244+1.307*0.2+(-3.6956)*x2+0.3486*x2.*x2; bb=polyfit(x2,yy5,2); xx5=5.25:0.05:7.25; yy51=polyval(bb,xx5); plot(xx5,yy51); grid on; figure(4)
yy6=29.1133+11.1342*0.2+(-7.608*x2)+0.6712*x2.*x2+(-1.4777)*x2*0.2; bb=polyfit(x2,yy6,2); xx6=5.25:0.05:7.25; yy61=polyval(bb,xx6); plot(xx6,yy61); grid on;
图1-7
图1-8
6. 模型的进一步改进
完全二次多项式模型：与1x 和2x 的完全二次多项式模型
02211223124152y x x x x x x ββββββε=++++++（6）
相比，模型（5）只少2
1x 项，我们不妨增加这一项，建立模型（10）。

这样做的好处之一是 MATLAB 统计工具箱有直接的命令rstool 求解，并且以交互式画面给出y 的估计值y ∧
和预 测空间。

代码为：
x=[x1' x2'];
rstool(x,y','quadratic')
结果为图1-9
图1-9
点击Export ，可以得到模型（6）的回归系数估计值为
β∧ =( 012345ββββββ∧∧∧∧∧∧
，
，，，，)=(32.0984,14.7436,-8.6367,-2.1038,1.1074,0.7594) 所以回归模型为：
Y=32.0984+14.7436*x1-8.6367*x2-2.1038*x1*x2+1.10742
1x +0.75942
2x
10.2软件开发人员的薪金
一家技术公司人事部门欲建立模型研究薪金与资历、管理责任、教育程度的关系，分析
人事策略的合理性，作为新聘用人员薪金的参考。

为此，研究人员收集了46名软件开发人员的档案资料，如表2-1，其中资历一列指从事专业工作的年数，管理一列中1表示管理人员，0表示非管理人员，教育一列中1表示中学程度，2表示大学程度，3表示更高程度（研究生）
一、问题重述
研究人员收集了46名软件开发人员的档案资料，以这资料建立模型研究薪金与资历、管理责任、教育程度的关系，分析人事策略的合理性，作为新聘用人员薪金的参考
二、问题分析
按照常识，薪金自然随着资历（年）的增长而增加，管理人员的薪金应高于非管理人员，教育程度越高薪金也越高
三、模型假设
1.建立薪金与资历，管理责任，教育程度之间的多元线性回归模型
2.利用matlab的统计工具箱计算回归系数及置信区间
3.在上述模型中增加管理与教育的交互项，建立新的回归模型
4.利用matlab的统计工具箱计算回归系数及置信区间并与上面结果比较得出结论
四、符号约定
对于问题，在符合题意并且与实际情况较吻合的情况下，薪金记作y，资历（年）记作x1,为了表示是
否非管理人员，定义x2={1,管理人员
0, 其它
，为了表示3种教育程度，定义x3={
1,中学
0,其它
，x4={
1,大学
0,其它
，
这样，中学用x3=1，x4=0来表示，大学用x3=0，x4=1表示，研究生则用x3=0,x4=0表示。

五、模型的建立与求解
1.基本模型
根据假设，薪金y与资历x1，管理责任x2，教育程度x3，x4之间的多元线性回归方程为：
y=a0+a1x1+a2x2+a3x3+a4x4+ （1）
2.模型求解
直接利用matlab统计工具箱中的命令regress求解回归系数估计值及其置信区间（置信水平α=0.05）、检验统计量2
R,F，p，2s，代码为：
x1=[111112222333344445556666788881010101011111212131314151616161720 ]';
x2=[1010010000111010000101011011000111001010110000]';
x3=[1000000100100100000110001010110000100100000101]';
x4=[0001011001010000110000110100001011010011011010]';
y=[13876116081870111283117672087211772105351219512313149752137119800114172026313 23112884132451367715965123662135213839228841697814803174042218413548144671594223174 2378025410148611688224170159902633017949256852783718838174831920719346]';
x0=ones(46,1);x=[x0x1x2x3x4];
[b,bint,r,rint,stats]=regress(y,x,0.05);
x0=ones(46,1);x=[x0x1x2x3x4];
[b,bint,r,rint,stats]=regress(y,x,0.05);
b=vpa(b,8)
bint=vpa(bint,8)
stats=vpa(stats,8)
结果如图2-1
图(2-1)
即模型（1
3.结果分析
从表2-2知2R=0.975，即因变量（薪金）的95.7%可由模型确定，F值远远超过F的检验的临界值，p
远小于α，因而模型（1）从整体来看是可用的。

比如，利用模型可以估计（或预测）一个大学毕业，有2年资历，费管理人员的薪金为：
y1=a0+a1*x1+a2*x2+a3*x3+a4*x4+ε=12272
模型中各个回归系数的含义可初步解释如下：x1的系数为546，说明资历增加1年薪金增长546；x2的系数为6883，说明管理人员薪金多6883；x3的系数为-2994，说明中学程度薪金比更高的少2994；x4的系数为148，说明大学程度薪金比更高的多148，但是应该注意到a4置信区间包含零点，说明这个系数的解释不可靠的。

需要指出，以上解释是就平均值来说，并且，一个因素改变引起的因变量的变化量，都是在其他因素不变的条件下成立的。

4.进一步的讨论
a4的置信区间包含零点，说明基本模型（1）存在缺点。

为了寻找改进的方向，常用残差分析方法（残差ε指薪金的实际值y与用模型估计的薪金y1之差，是模型（1）中随机误差ε的估计值，这里用了一个符号）。

我们将影响因素分成资历与管理-教育组合两类，管理-教育组合的定义如表三：
表2-3 管理-教育组合
为了对残差进行分析，下面用matlab绘图，代码为：
x5=[2 5 6 3 5 4 3 1 5 3 2 4 6 1 6 5 3 3 5 2 1 6 3 4 2 3 2 6 1 1 3 6 4 4 1 3 6 1 4 3 6 4 3 1 3 1]';
[b,bint,r,rint,stats]=regress(y,x,0.05);
subplot(2,2,1);
plot(x1,r,'r+');title('模型（1）随机误差与x1的关系');
subplot(2,2,2);
plot(x5,r,'b+');title('模型（1）随机误差与x2-x3，x4组合间的关系');
结果如图2-2
图（2-2）
从图1看，残差大概分成3个水平，这是由于6种管理—教育组合混合在一起，在模型中未被正确反映的结果，、；从图2看，对于前4个管理—教育组合，残差或者全为正，或者全为负，也表明管理—教育组合在模型中处理不当。

在模型（1）中国管理责任和教育程度是分别起作用的，事实上，二者可能起着交互作用，如大学程度的管理人员的薪金会比二者分别的薪金制和高一点。

以上分析提醒我们，应在基本模型（1）中增加管理x2与教育x3，x4的交互项，建立新的回归模型。

5.更好的模型
增加x2与x3，x4的交互项后，模型记作
y=a0+a1*x1+a2*x2+a3*x3+a4*x4+a5*x2*x3+a6*x3*x2+ （2）输入代码：
x0=ones(46,1);x=[x0x1x2x3x4x2.*x3x2.*x4];
[b,bint,r,rint,stats]=regress(y,x,0.05);b=vpa(b,8)
bint=vpa(bint,8)
stats=vpa(stats,8)
运行结果如图2-3
图（2-3）
即模型（2
表2-4 模型（2）计算结果
由表四可知，模型（2）的
2
R和F值都比模型（1）有所改进，并且所有回归系数的置信区间都不含
零点，表明模型（2）是完全可用的。

再与模型（1）类似，做残差分析图，程序为：
subplot(2,2,1);
plot(x1,r,'r+');
title('模型（2）随机误差与x1的关系');
subplot(2,2,2);
plot(x5,r,'b+');
title('模型（2）随机误差与x2-x3，x4组合间的关系');
结果为图2-4
图（2-4）
如图可以看出，已经消除了图2-2中的不正常现象，这也说明了模型（2）的适用性。

6.去掉异常数据，进一步优化模型：
图2-4中还可以发现一些异常点：具有10年资历，大学程度的管理人员（从表2-1中可以查出是33号），他的实际薪资明显低于模型的预估值，也明显低于与他有类似经历的其他人的薪金，这
可能是由于我们未知的原因造成的，为了使个别的数据不致影响整个模型，应将这个异常数据去
表2-5 模型（2）去掉异常数据后的计算结果
残差分析图代码为：
x1=[1 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 6 6 6 6 7 8 8 8 8 10 10 10 11 11 12 12 13 13 14 15 16 16
16 17 20]';
x2=[1 0 1 0 0 1 0 0 0 0 1 1 1 0 1 0 0 0 0 1 0 1 0 1 1 0 1 1 0 0 0 1 1 0 0 1 0 1 0 1 1 0 0 0 0]';
x3=[1 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 1 1 0 0 0 1 0 1 0 1 1 0 0 0 1 0 0 1 0 0 0 0 0 1 0 1]';
x4=[0 0 0 1 0 1 1 0 0 1 0 1 0 0 0 0 1 1 0 0 0 0 1 1 0 1 0 0 0 0 1 0 1 0 1 0 0 1 1 0 1 1 0 1 0]';
x5=[2 5 6 3 5 4 3 1 5 3 2 4 6 1 6 5 3 3 5 2 1 6 3 4 2 3 2 6 1 1 3 6 4 1 3 6 1 4 3 6 4 3 1 3 1]';
y=[13876 11608 18701 11283 11767 20872 11772 10535 12195 12313 14975 21371 19800 11417 20263 13231 12884 13245 13677 15965 12366 21352 13839 22884 16978 14803 17404 22184 13548 14467 15942 23174 25410 14861 16882 24170 15990 26330 17949 25685 27837 18838 17483 19207 19346]'; x0=ones(45,1);
x=[x0 x1 x2 x3 x4 x2.*x3 x2.*x4];
[b,bint,r,rint,stats]=regress(y,x,0.05);
vpa(b,8);vpa(bint,8);vpa(stats,8);
;plot(x1,r,'r+');title('模型（2）去异常数据后随机误差与x1间的关系');
得到的结果如图2-5
plot(x5,r,'b+');title('模型（2）去异常数据后随机误差与x2-x3，x4组合间的关系');
得到的结果如图2-6
7.模型应用
对于回归模型（2），用去掉模型异常数据（33号）后估计出的系数，得到的结果是满意的。

作为这个模型的应用之一，不妨用它来"制订”6种管理—教育组合人员的“基础”薪金（即资历为零的薪金，当然，这也是平均意义上的）。

利用模型（2）和表2-5容易得到表2-6：
图（2-5）
图（2-6）
可以看出，大学程度的管理人员的薪金比研究生程度的管理人员的薪金高，而大学程度的非管理人员的薪金比研究生程度的管理人员的薪金高，而大学程度的非管理人员的薪金比研究生的非管理人员的薪金略低。

当然，这是根据这家公司实际数据建立的模型得到的结果，并不具普遍性
实验小结
这是本学期数学建模实验课最后一次实验报告了，经过半年多的学习，我初步掌握了MATLAB这个软件的使用方法。

也大致了解了这个软件的功能，知道如何利用软件去解决一些问题。

但我深知我所掌握的东西只是这个软件的九牛一毛而已，虽然课程即将结束，但这个软件和使用技巧不应该被我束之高阁，而应该在平时自己多琢磨，多钻研，熟练地使用更多的功能与技巧，解决更复杂和更重要的问题。