数学建模举例

10.1牙膏的销售量某大型牙膏制造企业为了更好地拓展产品市场，有效地管理库存，公司董事会要求销售部门根据市场调查，找出公司生产的牙膏销售量与销售价格、广告投入等之间的关系，从而预测出在不同价格和广告费用下的销售量。为此，销售部的研究人员收集了过去30个销售周期（每个销售周期为4周）公司生产的牙膏的销售量、销售价格、投入的广告费用，以及同期其它厂家生产的同类牙膏的市场平均销售价格，见表1-1（其中价格差指其它厂家平均价格与公司销售价格之差）。试根据这些数据建立一个数学模型，分析牙膏销售量与其它因素的关系，为制订价格策略和广告投入策略提供数量依据表1-1牙膏销售量与销售价格、广告费用等数据销售周期公司销售价格(元)其他厂家平均价格(元)价格差(元)广告费用(百万元)销售量(百万支)13.853.80-0.055.57.3823.754.000.256.758.5133.704.300.607.259.5243.603.700.005.507.5053.603.850.257.009.3363.63.800.206.508.2873.63.750.156.758.7583.83.850.055.257.8793.83.65-0.155.257.10103.854.000.156.008.00113.904.100.206.507.89123.904.000.106.258.15133.704.100.407.009.10143.754.200.456.908.86153.754.100.356.808.90163.804.100.306.808.87173.704.200.507.109.26183.804.300.507.009.00193.704.100.406.808.75203.803.75-0.056.507.95213.803.75-0.056.257.65223.753.65-0.106.007.27233.703.900.206.508.00243.553.650.107.008.50253.604.100.506.808.75263.704.250.606.809.21273.753.65-0.056.508.27283.753.750.005.757.67293.803.850.055.807.93303.704.250.556.809.26一、 问题重述根据过去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的值有比较明显的线性增长趋势，图中的直线是用线性模型：（1）拟合的（其中是随机误差）。而在图2中，当x2增大时，y有向上弯曲增加的趋势，图中的曲线是用二次函数模型:（2）综合上面的分析，结合模型（1）和（2）建立如下的回归模型:（3）（3）式右端的x1和x2称为回归变量（自变量），是给定价格差x1，广告费用x2时，牙膏销售量y的平均值，其中的参数称为回归系数，由表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）、检验统计量,F，p，得结果见表1-2，参数参数估计值参数置信区间17.32445.7282,28.92061.30700.6829,1.9311-3.6956-7.4989,0.10770.34860.0379,0.6594=0.9054 F=82.9409 p0.0001 =0.0490表1-2 模型（3）计算结果图（1-3）3. 结果分析表1-2显示，=0.9054指因变量y（销售量）的90.54%可由模型决定，F值远远超过F检验的临界值，P远小于，因而模型（3）整体来看是可用的表1-2的回归系数中的置信区间包含零点，表示回归变量（对因变量y的影响）是不太显著的，但由于是显著的，我们仍将留在模型中4. 销售量预测将回归系数的估计值带入模型（3），即可预测公司未来某个销售周期牙膏的销售量y，预测值记为，得到模型（3）的预测方程：（4）只需要知道该销售周期的价格差x1和投入的广告费用x2,就可以计算预测值。5. 模型改进模型（3）中回归变量x1和x2对因变量y的影响是相互独立的，即牙膏销售量y的均值与广告费用x2的二次关系由回归系数和确定，而不依赖于价格差x1，同样的，y的均值与x1的线性关系由回归系数确定，而不依赖于x2。根据直觉和经验可以猜想，x1和x2之间的交互作用会对y有影响，不妨简单地用x1，x2的乘积代表它们的相互作用，于是将模型（3）增加一项，得到:（5）在这个模型中，y的均值与2x的二次关系为，由系数确定，并依赖于价格差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）计算结果即为表1-3参数参数估计值参数置信区间29.113313.7013,44.525211.13421.9778,20.2906-7.6080-12.6932,-2.52280.67120.2538,1.0887-1.4777-2.8518,-0.1037=0.9209 F=72.7771 P0.0001 =0.0426表1-3 模型（5）计算结果表3与表2的结果相比，有所提高，说明模型（5）比模型（3）有所进步。并且，所有参数的置信区间，特别是X1，X2的交互作用项X1X2的系数的置信区间不包含零点，所以有理由相信模型（5）比模型（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 onfigure(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）中牙膏销售量的均值与广告费用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-86. 模型的进一步改进完全二次多项式模型：与1x和2x的完全二次多项式模型（6）相比，模型（5）只少项，我们不妨增加这一项，建立模型（10）。这样做的好处之一是MATLAB统计工具箱有直接的命令rstool求解，并且以交互式画面给出y的估计值和预测空间。代码为：x=x1 x2;rstool(x,y,quadratic)结果为图1-9图1-9点击Export，可以得到模型（6）的回归系数估计值为 =( )=(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.1074+0.759410.2软件开发人员的薪金一家技术公司人事部门欲建立模型研究薪金与资历、管理责任、教育程度的关系，分析人事策略的合理性，作为新聘用人员薪金的参考。为此，研究人员收集了46名软件开发人员的档案资料，如表2-1，其中资历一列指从事专业工作的年数，管理一列中1表示管理人员，0表示非管理人员，教育一列中1表示中学程度，2表示大学程度，3表示更高程度（研究生）表2-1 软件开发人员的薪金与资历、管理责任、教育程度的关系编号薪金资历管理教育011387611102116081030318701113041128310205117671030620872212071177220208105352010912195203101231330211149753111221371312131980031314114174011520263413161323140317128844021813245502191367750320159655112112366601222135261323138396022422884612251697871126148038022717404811282218481329135488013014467100131159421002322317410133323780101234254101112351486111013616882120237241701213381599013013926330131240179491402412568515134227837161243188381602441748316014519207170246193642001一、 问题重述研究人员收集了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）、检验统计量,F，p，代码为：x1=111112222333344445556666788881010101011111212131314151616161720;x2=1010010000111010000101011011000111001010110000;x3=1000000100100100000110001010110000100100000101;x4=0001011001010000110000110100001011010011011010;y=13876116081870111283117672087211772105351219512313149752137119800114172026313231128841324513677159651236621352138392288416978148031740422184135481446715942231742378025410148611688224170159902633017949256852783718838174831920719346;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）的计算结果是表2-2参数参数估计值参数置信区间a01103210258,11807a1546484,608a268836248,7517a3-2994-3826,-2162a4148-636,931=0.957 F=226 p0.0001 =1.057*106表2-2 模型（1）计算结果3. 结果分析从表2-2知=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）中随机误差的估计值，这里用了一个符号）。我们将影响因素分成资历与管理-教育组合两类，管理-教育组合的定义如表三：组合123456管理010101教育112233表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参数参数估计值参数置信区间a01120411044,11363a1497486,508a270486841,7255a3-1727-1939,7255a4-348-545,-152a5-3071-3372,-2769a618361571,2101=0.9988 F=5545 P0.0001 =3.0047*104表2-4 模型（2）计算结果由表四可知，模型（2）的和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）重新估计回归系数得到结果如表2-5参数参数估计值参数置信区间a01120011139,11261a1498494，503a270416962,7120a3-1737-1818，-1656a4-356-431，-281a5-3056-3171，-2942a619971894,2100=0.9988 F=36701P0.0001 =4.347*103表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