福建农林大学计算机与信息学院
（数学类课程）
课程名称：姓名：系：专业：年级：学号：指导教师：职称：实验报告
数学模型 苏志东 数学 数学与应用数学
2014级
姜永 副教授
2016年 6月12日
实验项目列表
序号 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 实验项目名称 数学规划模型建立及其软件求解 数据插值与数据拟合应用 统计回归模型及其软件求解 成绩 指导教师 姜 永 姜 永 姜 永
福建农林大学计算机与信息学院数学类实验报告（一）
系： 数学 专业： 数学与应用数学 年级： 2014级 姓名： 学号： 3 实验课程： 数学模型 实验室号： 明南附203 实验设备号：实验时间： 2016/6/6 指导教师签字：成绩：
1．实验项目名称：数学规划模型建立及其软件求解 2．实验目的和要求：
了解数学规划的的基本理论和方法，并用于建立实际问题的数学规划模型；会用LINGO软件解数学规划问题并对结果加以分析应用。
3．实验使用的主要仪器设备和软件：联想启天M430E电脑；LINGO12.0或以上版本。 4．实验的基本理论和方法：
一般地，数学规划模型可表述成如下形式：
Minz?f(x)
xs.t.gi(x)?0,i?1,2,...,m.
g(x)?0(i?1,2,...,m)为约束条件。
其中f(x)表示目标函数，i LINGO用于解决二次规划、线性规划以及非线性规划问题，同时可以求解线性或非线性方程（组）。LINGO的最大特色在于通过高运行速度解决优化模型中的决策变量的整数取值问题。
线性优化求解程序通常使用单纯性算法，可以使用LINGO的内点算法解决大规模规划问题。非线性规划可通过迭代求解一系列线性规划求解。
5．实验内容与步骤：
问题一：某公司将3种不同含硫量的液体原料（分别记为甲、乙、丙）混合生产两种产品（分别记为A，B），按照生产工艺的要求，原料甲、乙必须首先倒入混合池中混合，混合后的液体再分别与原料丙混合生产A，B．已知原料甲，乙，丙的含硫量分别是3%，1%，2%，进货价格分别为6千元/ t，16千元/ t，10千元/t，产品A，B的含硫量分别不能超过2.5%，1.5%，售价分别为9千元/t，15千元/t，根据市场信息，原料甲、乙、丙的供应量都不能超过500t；产品A，B的最大市场需求量分别为100t ,200t．
(1) 应如何安排生产？
(2) 如果产品A的最大市场需求量增长为600t，应如何安排生产？
(3) 如果乙的进货价格下降为12千元/t，应如何安排生产？分别、对（1）、（2）两种情况进行讨论．
解答：
(1)问题分析
根据题目要求，不难想到，这个问题的目标是使公司获利最大，要做的决策就是生产计划，即生产多少产品A和产品B ，限制条件有：原料供应、市场需求、不同含硫量生产不同的产品。根据这些条件，利用lingo软件，求出最终决策。 基本模型
决策变量：设用（i=甲，乙，丙；j=A,B）表示用第i种原料用于生产产品j，将i=甲,乙,丙转换为i=1,2,3,j=A,B转换为j=1,2.
目标函数:设公司获利为z千元,则有:
3
z?9?xi1?15?xi2?6?x1j?16?x2j?10?x3ji?1i?1j?1j?1j?133222
约束条件
原料供应:原料i(i=1,2,3)均不超过500t,则
?xj?12ij?500 (i=1,2,3)
市场需求:产品A、B的需求量分别为100t、200t，则有：
?3xi1?100???i?1?3??xi2?200??i?1
含硫量：根据甲乙混合比例，有：x11:x21?x12:x22，由生产不同产品含硫量百分比，有：
3%x11?1%x21?2%x31?1.5%??2.5%?x?x?x?112131?3%x12?1%x22?2%x32?1%??1.5%?x?x?x122232?
终上所述，有：maxz?933222
?xi?1i1?15?xi2?6?x1j?16?x2j?10?x3j
i?1j?1j?1j?1?xj?12ij?500 (i=1,2,3)
?x?xi?1i?133i1?100?200
i21.5%?3%x11?1%x21?2%x31?2.5%x11?x21?x313%x12?1%x22?2%x321%??1.5%x12?x22?x32
对上述式子进行调整，并利用lingo软件，可求解出最优解。
Lingo程序为：
max=9*(x11+x21+x31)+15*(x12+x22+x32)-6*(x11+x12)-16*(x21+x22)-10*(x31+x32); 0.5*x11-1.5*x21-0.5*x31<=0; 1.5*x11-0.5*x21+0.5*x31>0; 1.5*x12-0.5*x22+0.5*x32<=0; 2*x12+x32>0;
4
x11*x22-x21*x12=0; x11+x12<=500; x21+x22<=500; x31+x32<=500;
x11+x21+x31<=100; x12+x22+x32<=200;
程序运行结果如下：
Objective value: 400.0000
Variable Value X11 0.000000 X21 0.000000 X31 0.000000 X12 0.000000 X22 100.0000 X32 100.0000
结果分析：根据结果显示，最优解为用100t的乙原料和100t的丙原料混合，生成200t产品B，所以目标函数最优解为40万元（400千元）。
（2）本小题的解法与（1）基本一致，只需要将约束条件改变为，相应的代码由x11+x21+x31<=100改为x11+x21+x31<=600，并代入程序计算，便可求解出结果。 程序运行结果如下：
Objective value: 600.0000
Variable Value X11 300.0000 X21 0.000000 X31 300.0000 X12 0.000000 X22 0.000000 X32 0.000000
结果分析：根据结果显示，最优解为用300t的甲原料和300t的丙原料混合，生成600t产品A所以目标函数最优解为60万元（600千元）。
（3）将乙的进货价格下降为12千元/t，只需修改一下目标函数值和约束条件即可。针对问题（1）来说，只需将目标函数z?933?xi?123i1?15?xi2?6?x1j?16?x2j?10?x3j改为
i?12j?1j?1j?123222对应的程序修改一下，即可得到新的求解结果。 z?9?xi1?15?xi2?6?x1j?12?x2j?10?x3j，
i?1i?1j?1j?1j?1程序运行结果如下：
Objective value: 900.0000 Variable Value Reduced Cost
X11 0.000000 0.000000 X21 0.000000 0.000000 X31 0.000000 0.000000 X12 50.00000 0.000000
5
X22 150.0000 0.000000 X32 0.000000 1.000000
结果分析：根据结果显示，最优解为用50t的甲原料和150t的乙原料混合，生成200t产品B，所以目标函数最优解为90万元（900千元）。
问题二：某造船厂需要决定下四个季度的帆船生产量。下四个季度的帆船需求量分别是40条、60条、75条和25条，这些需求必须按时满足。每个季度正常的生产能力是40条帆船，每条船的生产费用为40万元。如果加班生产，每条船的生产费用为45万元。每个季度末，每条船的库存为2万元。假定生产提前期为0，初始库存为10条船。如何安排生产可使总费用最小？(LINGO程序要求利用集合语言编写)
解答： 建立模型
设四个季度轮船的需求量分别为
DEM(I),I?1,2,3,4;
四个季度正常生产的产量分别为
RP(I),I?1,2,3,4；
四个季度加班生产的产量分别为
OP(I),I?1,2,3,4 ；
四个季度轮船的总量分别为
ALL(I),I?1,2,3,4
根据题意和约束条件可以建立以下模型： 目标函数：
?(40*RP(I)?45*OP(I)?2*(ALL(I)?DEM(I)))
I?14约束条件由题意依次为
1、每季度正常生产能力是40条船，即
I?1,2,3,4,应有RP(I)??40;
2、需求量限制:
I?1,2,3,4,应有ALL(I)??DEM(I);
模型求解
利用题目所给数据，将所建立的目标函数以及限制条件输入LINGO:
模型代码如下： sets:
SIJI/1..4/:DEM,RP,OP,ALL; endsets data:
DEM=40 60 75 25; enddata
6
ALL(1)=10+RP(1)+OP(1);
ALL(2)=ALL(1)-DEM(1)+RP(2)+OP(2); ALL(3)=ALL(2)-DEM(2)+RP(3)+OP(3); ALL(4)=ALL(3)-DEM(3)+RP(4)+OP(4);
min=@sum(SIJI(I):40*RP(I)+45*OP(I)+2*(ALL(I)-DEM(I))); @for(SIJI(I):RP(I)<=40);
@for(SIJI(I): ALL(I)>=DEM(I)); end
点击运行按钮得试验结果如下： Global optimal solution found.
Objective value: 7845.000
Variable Value Reduced Cost DEM( 1) 40.00000 0.000000 DEM( 2) 60.00000 0.000000 DEM( 3) 75.00000 0.000000 DEM( 4) 25.00000 0.000000 RP( 1) 40.00000 0.000000 RP( 2) 40.00000 0.000000 RP( 3) 40.00000 0.000000 RP( 4) 25.00000 0.000000 OP( 1) 0.000000 2.000000 OP( 2) 10.00000 0.000000 OP( 3) 35.00000 0.000000 OP( 4) 0.000000 5.000000 ALL( 1) 50.00000 0.000000 ALL( 2) 60.00000 0.000000 ALL( 3) 75.00000 0.000000 ALL( 4) 25.00000 0.000000 结果分析：RP(1)?40,RP(2)?40,RP(3)?40,RP(4)?25;
OP(1)?0,OP(2)?10,OP(3)?35,OP(4)?0。
所以须这样安排生产：第一季度需生产40条,无需加班；第二季度需生产出50条,其中有10条是加班生产的；第三季度需生产出75条,其中35条是加班生产的；第四季度需生产出25条，无需加班；最小总费用为7845万元。
问题三：某人事部欲安排四个人到四个不同的岗位工作，每个岗位一个人，经考核四人在不同岗位的成绩如下表，应如何安排他们的工作才能使总成绩最好？(LINGO程序要求利用集合语言编写)
工作 A B C D 人员 甲 乙 85 95 91 88 70 78 90 93 7
丙 丁 82 86 84 89 79 81 90 88
解答：记甲乙丙丁分别为人员i=1,2,3,4；记工作A、B、C、D分别为j=1,2,3,4.记人员i的
第j种工作的最好成绩为基本模型
。
min z =
约束条件:
,i=1,2,3,4 ,j=1,2,3,4
对上述式子进行调整，并利用lingo软件，可求解出最优解。 Lingo程序为：
model: sets:
person/1..4/; position/1..4/;
link(person,position):c,x; endsets data:
c=85,91,70,90, 95,88,78,93, 82,84,79,90, 86,89,81,88; enddata
max=@sum(link:c*x);
@for(person(i):@sum(position(j):x(i,j))<=1); @for(position(i):@sum(person(j):x(j,i))=1); @for(link:@bin(x)); end
程序运行结果如下：
Global optimal solution found.
Objective value: 357.0000 Objective bound: 357.0000 Infeasibilities: 0.000000 Extended solver steps: 0 Total solver iterations: 0
Model Class: PILP
8
Total variables: 16 Nonlinear variables: 0 Integer variables: 16
Total constraints: 9 Nonlinear constraints: 0
Total nonzeros: 48 Nonlinear nonzeros: 0
Variable Value Reduced Cost C( 1, 1) 85.00000 0.000000 C( 1, 2) 91.00000 0.000000 C( 1, 3) 70.00000 0.000000 C( 1, 4) 90.00000 0.000000 C( 2, 1) 95.00000 0.000000 C( 2, 2) 88.00000 0.000000 C( 2, 3) 78.00000 0.000000 C( 2, 4) 93.00000 0.000000 C( 3, 1) 82.00000 0.000000 C( 3, 2) 84.00000 0.000000 C( 3, 3) 79.00000 0.000000 C( 3, 4) 90.00000 0.000000 C( 4, 1) 86.00000 0.000000 C( 4, 2) 89.00000 0.000000 C( 4, 3) 81.00000 0.000000 C( 4, 4) 88.00000 0.000000 X( 1, 1) 0.000000 -85.00000 X( 1, 2) 1.000000 -91.00000 X( 1, 3) 0.000000 -70.00000 X( 1, 4) 0.000000 -90.00000 X( 2, 1) 1.000000 -95.00000 X( 2, 2) 0.000000 -88.00000 X( 2, 3) 0.000000 -78.00000 X( 2, 4) 0.000000 -93.00000 X( 3, 1) 0.000000 -82.00000 X( 3, 2) 0.000000 -84.00000 X( 3, 3) 0.000000 -79.00000 X( 3, 4) 1.000000 -90.00000 X( 4, 1) 0.000000 -86.00000 X( 4, 2) 0.000000 -89.00000 X( 4, 3) 1.000000 -81.00000 X( 4, 4) 0.000000 -88.00000
Row Slack or Surplus Dual Price
9
1 357.0000 1.000000 2 0.000000 0.000000 3 0.000000 0.000000 4 0.000000 0.000000 5 0.000000 0.000000 6 0.000000 0.000000 7 0.000000 0.000000 8 0.000000 0.000000 9 0.000000 0.000000
结果分析：让甲到B岗位工作，乙到A岗位工作，丙到D岗位工作，丁到C岗位工作可以使总成绩最好，为357。
6．实验心得（质疑、建议）：
10
福建农林大学计算机与信息学院数学类实验报告（二）
系： 数学 专业： 数学与应用数学 年级： 2014级 姓名： 学号： 3 实验课程： 数学模型 实验室号： 明南附203 实验设备号：实验时间： 2016/6/6 指导教师签字：成绩：
1．实验项目名称：数据插值与数据拟合应用 2．实验目的和要求：
理解数据插值与数据拟合的理论和方法，会使用MATLAB进行数据插值与数据拟合，能够使用MATLAB解决一些关于数据插值与数据拟合的应用问题。
3．实验使用的主要仪器设备和软件：联想启天M430E电脑；MATLAB2010或以上版本。 4．实验的基本理论和方法：
4.1插值与拟合
在实际工程应用和科学实际和科学实践中，经常需要寻求两个（或多个）变量间的关系，而实际却只能通过观测得到一些离散的数据点。针对分散的数据点，运用某种数学方法确定两个（或多个）变量间的函数关系，这个过程称为数据插值或数据拟合。
假设x为自变量，y为因变量，函数关系为y?f(x)（待定）。现给定一组点（xi,yi）,(i?1,2,...,n)，然后构造一个简单函数P（x）作为函数y?f(x)的近似表达式，即
y?f(x)?P(x)（1）
对式（1），若满足
P（xi）?f(xi)?yi,i?1,2,...,n（2）
这类问题称为插值问题。
式（2）要求所求的函数曲线通过已知的数据点，若不要求P（x）通过所有数据点（xi,yi）,(i?1,2,...,n)，而是要求曲线在某种准则下整体与所给的数据点尽量接近，如按最小二乘法要求??达到最小，而得到P（x），此类问题称为拟合问题。 4.2最小二乘法
给定平面上一组点(xi,yi)(i=1,2,3,...,n),作曲线拟合有多种方法,其中最小二乘法是常用的一种。最小二乘法的原理是：求f(x),使??
?(f(x)?y)iii?1n2?(f(x)?y)达到最小。拟合时选取一定的拟合函数形式。
2iii?1n11
5．实验内容与步骤：
问题一（插插值问题）：有一组数据如下，试用不同的插值方法分别计算x1?1.56，x2?6.23所对应的近似值。
x y
解答： （1）线性插值 x=1:6;
y=[1.0000 1.2599 1.4422 1.5874 1.7100 1.8171]; xi=1:0.1:6;
yi=interp1(x,y,xi,'linear'); plot(xi,yi,'k',x,y,'o') axis tight x0=1.56;
y0=interp1(x,y,x0,'linear') y0 = 1.1455 x0=6.23;
y0=interp1(x,y,x0,'linear') y0 = NaN
1.81.71.61.51.41.31.21.1111 2 3 4 5 6
1.0000 1.2599 1.4422 1.5874 1.7100 1.8171
1.522.533.544.555.56
12
（2）最近邻点插值 x=1:6;
y=[1.0000 1.2599 1.4422 1.5874 1.7100 1.8171]; xi=1:0.1:6;
yi=interp1(x,y,xi,'nearest'); plot(xi,yi,'k',x,y,'o') axis tight x0=1.56;
y0=interp1(x,y,x0,'nearest') y0 = 1.2599 x0=6.23;
y0=interp1(x,y,x0,'nearest') y0 = NaN
1.81.71.61.51.41.31.21.1111.522.533.544.555.56
（3）三次样条函数插值 x=1:6;
y=[1.0000 1.2599 1.4422 1.5874 1.7100 1.8171]; xi=1:0.1:6;
yi=interp1(x,y,xi,'spline'); plot(xi,yi,'k',x,y,'o')
13
axis tight x0=1.56;
y0=interp1(x,y,x0,'spline') y0 = 1.1579 x0=6.23;
y0=interp1(x,y,x0,'spline') y0 = 1.8403
1.81.71.61.51.41.31.21.1111.522.533.544.555.56
（4）三次函数插值 x=1:6;
y=[1.0000 1.2599 1.4422 1.5874 1.7100 1.8171]; xi=1:0.1:6;
yi=interp1(x,y,xi,'cubic'); plot(xi,yi,'k',x,y,'o') axis tight x0=1.56;
y0=interp1(x,y,x0,'cubic') y0 = 1.1560 x0=6.23;
14
y0=interp1(x,y,x0,'cubic') y0 = 1.8395
1.81.71.61.51.41.31.21.1111.522.533.544.555.56
问题二（给药问题）：一种新药用于临床之前，必须设计给药方案，即每次注射计量多大，间隔时间多长。
药物进入机体后随血液输送到全身，在这个过程中不断被吸收、分布、代谢、最终排除体外。药物在血液中的浓度，即单位体积血液中的药物含量，称为血药浓度。在最简单的一室模型中，将整个机体看作一个房室，称为中心室，室内的血药浓度是均匀的。快速静脉注射后，浓度立即上升，然后逐渐下降。当浓度太低时，达不到预期的治疗效果；当浓度太高时，又可能导致药物中毒或副作用太强。根据临床经验要求：每种药物有一个最小有效浓度c1和一个最大浓度c2。设计给药方案时，要使血药浓度保持在c1~c2之间，本问题设c1?10,c2?25(mg/ml)，而且本问题可视为一室模型。
设对某人用快速静脉注射方式—次注入某药物300mg后，在一定时刻t(h)采取血药，测得血药浓度c(mg/ml)如下表：
t c
0.25 19.21
0.5 18.15
1 15.36
1.5 14.10
2 12.89
3 9.32
4 7.45
6 5.24
8 3.01
请根据上述数据，利用房室模型和数据拟合方法确定给药方案。
解答：
（1）根据题目提供的数据及提示，为了更好地解决问题。我们可以假设：整个过程中血液容积不变，建立如下模型：
c(t)?依题意可知：初值
D0?kte V15
D0?300(mg),c(0)?20(mg/ml)，
从而我们可以知道
V?根据测得的浓度可知，在
时刻，
300?15(ml)20。
，由模型可得以下式子：
c(8)?解得：
300?8k?e 15k?0.24
根据以上计算，我们可以得到模型的初值：
（2）根据初值，通过MATLAB编程，具体如下： 首先建立M-文件，该文件命令为：curvefun1.m function f=curvefun1(x,tdata) f=x(1)/x(2)*exp(-x(3)*tdata);
然后在command windows(命令窗口)输入以下程序： tdata=[0.25 0.5 1 1.5 2 3 4 6 8];
cdata=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01]; x0=[300 15 0.24];
x=lsqcurvefit('curvefun1',x0,tdata,cdata) 根据程序得出以下结果：
（3）由已知可得：
c1?10(mg/ml),c2?25(mg/ml)
并根据题意可得： 首次注射量
每次注射量
间隔时间
??lnc2?lnc1ln25?ln10??3.7863k0.2420
6．实验心得（质疑、建议）：
福建农林大学计算机与信息学院数学类实验报告（三）
系： 数学 专业： 数学与应用数学 年级： 2014级
16
姓名： 学号： 3 实验课程： 数学模型
实验室号： 明南附203 实验设备号：实验时间： 2016/6/6 指导教师签字：成绩：
1．实验项目名称：统计回归模型及其软件求解 2．实验目的和要求：
了解回归分析的基本原理，掌握MATLAB实现的方法；学习应用回归模型解决实际问题。
3．实验使用的主要仪器设备和软件：联想启天M430E电脑；MATLAB2010或以上版本。 4．实验的基本理论和方法：
当人们对研究对象的内在特性和各因素间的关系有比较充分的认识时，一般用机理分析方法建立
数学模型。如果由于客观事物内部规律的复杂性及人们认识程度的限制，无法分析实际对象内在的因果关系，建立合乎机理规律的数学模型，那么通常的办法是搜集大量的数据，基于对数据的统计分析去建立模型，其中一类应用非常广泛的随机模型就是统计回归模型。
回归分析是研究一个变量Y与其他若干变量X之间相关关系的一种数学工具，它是在一组试验或观测数据的基础上，寻找被随机性掩盖了的变量之间的依存关系。粗略的讲，可以理解为用一种确定的函数关系去近似代替比较复杂的相关关系，这个函数称为回归函数，在实际问题中称为经验公式。回归分析所研究的主要问题就是如何利用变量X,Y的观察值，对回归函数进行统计推断，包括对它进行估计及检验与它有关的假设等。
Matlab命令：散点图：plot（x,y,’o’）回归工具箱：rstool
线性回归：[b,bint,r,rint,stats]=regress(y,x,alpha) 残差图：rcoplot(r,rint) 多项式回归：[p,S]=plotfit(x,y,m)
非线性回归：[beta,R,J]=nlinfit(x,y,’model’,bata0);nlparci(beta,R,J);nlintool 逐步回归：stepwise(x,y,inmodel,alpha)
5．实验内容与步骤：
问题一：下表列出了某城市18位35岁－44岁经理的年平均收入x1千元，风险偏好度x2和人寿保险额y千元的数据，其中风险偏好度是根据发给每个经理的问卷调查表综合评估得到的，它的数值越大，就越偏爱高风险。研究人员想研究此年龄段中的经理所投保的人寿保险额与年均收入及风险偏好度之间的关系。研究者预计，经理的年均收入和人寿保险额之间存在着二次关系，并有把握地认为风险偏好度对人寿保险额有线性效应，但对风险偏好度对人寿保险额是否有二次效应以及两个自变量是否对人寿保险额有交互效应，心中没底。
请你通过表中的数据通过试验来建立一个合适的回归模型，验证上面的看法，并给出进一步的分析。
17
序号 1 2 3 4 5 6 7 8 9 y 196 63 252 84 126 14 49 49 266 x1 66.290 40.964 72.996 45.010 57.204 26.852 38.122 35.840 75.796 x2 7 5 10 6 4 5 4 6 9 序号 10 11 12 13 14 15 16 17 18 y 49 105 98 77 14 56 245 133 133 x1 37.408 54.376 46.186 46.130 30.366 39.060 79.380 52.766 55.916 x2 5 2 7 4 3 5 1 8 6 解答：
基本模型1：由题目可知，经理的年均收入和人寿保险额之间存在着二次关系，且风险偏好度对人寿保险额有线性效应。综上所述，建立如下的回归模型
y??0??1x2??2x12 （1）
其中的参数?0,?1,?2是回归系数。
模型求解：利用MATLAB统计工具箱中的命令regress求解：
x1=[66.290 40.964 72.996 45.010 57.204 26.852 38.122 35.840 75.796 37.408 54.376 46.186 46.130 30.366 39.060 79.380 52.766 55.916]';
x12=x1.^2;
x2=[7 5 10 6 4 5 4 6 9 5 2 7 4 3 5 1 8 6]';
y=[196 63 252 84 126 14 49 49 266 49 105 98 77 14 56 245 133 133]'; x=[ones(18,1) x2 x12];
[b,bint,r,rint,stats]=regress(y,x); b,bint,stats,rcoplot(r,rint) 运行结果：
b =-41.8542 5.8123 0.0447 bint =-45.5187 -38.1897 5.2196 6.4050 0.0439 0.0455 stats =1.0e+003
0.0010 8.1851 0.0000 0.0066
18
Residual Case Order Plot8642Residuals0-2-4-6-8246810Case Number12141618 图1 模型1残差分布图
2结果分析：结果显示，R?1指因变量（人寿保险额）的100%可由模型确定，F值远远超过F检验的
临界值，p值为0，模型中的?0??41.8542,?1?5.8123,?2?0.0447，且它们的置信区间都不包含零点，因而模型从整体来看是可用的。但存在异常数据。
模型1改进：剔除第6组异常数据模型求解：
模型求解：利用MATLAB统计工具箱中的命令regress求解：
x1=[66.290 40.964 72.996 45.010 57.204 38.122 35.840 75.796 37.408 54.376 46.186 46.130 30.366 39.060 79.380 52.766 55.916]';
x12=x1.^2;
x2=[7 5 10 6 4 4 6 9 5 2 7 4 3 5 1 8 6]';
y=[196 63 252 84 126 49 49 266 49 105 98 77 14 56 245 133 133]'; x=[ones(17,1) x2 x12];
[b,bint,r,rint,stats]=regress(y,x); b,bint,stats,rcoplot(r,rint) 运行结果：
b =-40.8713 5.8314 0.0444 bint =-44.0484 -37.6942 5.3327 6.3302 0.0437 0.0451 stats = 1.0e+004 *
0.0001 1.0625 0.0000 0.0005
19
Residual Case Order Plot642Residuals0-2-4-6246810Case Number121416 图2 改进模型1的残差分布图 结果分析：改进后的模型的F有所提高,s有所下降。而且预测区间长度短一些，精度提高了。所以改进后的模型为y??40.8713?5.8314x2?0.0444x12，比之前的模型更优。
进一步讨论：根据直觉和经验猜想，经理的年收入可能会对人寿保险额有线性效应，风险偏好度可能会对人寿保险额会有二次效应，以及两个自变量可能会对人寿保险额有交互效应。于是，将模型改为
2 （2） y??0??1x1??2x2??3x1x2??4x12??5x22模型2求解：利用MATLAB统计工具箱中的命令regress求解：
x1=[66.290 40.964 72.996 45.010 57.204 26.852 38.122 35.840 75.796 37.408 54.376 46.186 46.130 30.366 39.060 79.380 52.766 55.916]';
x2=[7 5 10 6 4 5 4 6 9 5 2 7 4 3 5 1 8 6]'; x1x2=x1.*x2; x12=x1.^2; x22=x2.^2;
y=[196 63 252 84 126 14 49 49 266 49 105 98 77 14 56 245 133 133]'; x=[ones(18,1) x1 x2 x1x2 x12 x22]; [b,bint,r,rint,stats]=regress(y,x); b,bint,stats,rcoplot(r,rint) 运行结果：
b =-65.3856 1.0172 5.2171 -0.0196 0.0358 0.1662 bint =-78.7266 -52.0447 0.5202 1.5141 2.2785 8.1558 -0.0501 0.0109 0.0310 0.0406 -0.0956 0.4279 stats =1.0e+03 *
0.0010 7.1102 0.0000 0.0030
20
2模型2分析：?3,?5的置信区间包含零点，表明回归变量x1x2,x2对因变量y的影响不是太显著，但x1是显著的，所以将x1保留在模型中。 模型2改进：将模型改为：
y??0??1x1??2x2??3x12
模型2求解：利用MATLAB统计工具箱中的命令regress求解：
将x=[ones(18,1) x1 x2 x1x2 x12 x22];改为x=[ones(18,1) x1 x2 x12]; 运行结果：
b =-62.3489 0.8396 5.6846 0.0371 bint = -73.5027 -51.1952 0.3951 1.2840 5.2604 6.1089 0.0330 0.0412 stats =1.0e+04 *
0.0001 1.1070 0.0000 0.0003
Residual Case Order Plot642Residuals0-2-4-6246810Case Number12141618 图3 模型2的残差分布图
2模型2分析：结果显示，R?1指因变量（人寿保险额）的100%可由模型确定，F值远远超过F检验
的临界值，p值为0，模型中的?0??62.3489,?1?0.8396,?2?5.6846,?3?0.0371，且它们的置信
x1从整体来看是可用的。区间都不包含零点，因而模型y??62.3489?0.8396x1?5.6846x2?0.0371 x1。 综上，改进后的模型2为最满意的模型。即y??62.3489?0.8396x1?5.6846x2?0.0371^2^2
问题二：某一个医药公司的新药研究部门为了掌握一种新止痛剂的疗效，设计了一个药物实验，
给患者有同种病痛的病人使用这种新止痛剂的以下4种剂量中的某一个：2g，5g，7g，10g，并记录每个病人病痛明显减轻的时间（以分钟计）。为了解新药的疗效与病人性别和血压有什么关系，试验过
21
程中研究人员把病人按性别及血压的低、中、高三档平均分配来进行测试。通过比较每个病人血压的历史数据，从低到高分成3组，分别记作0.25，0.50，0.75。实验结束后，公司的记录结果见下表（性别以0表示女，1表示男）。
请你利用统计回归方法为公司建立一个模型，根据病人用药的剂量、性别和血压级别，预测出服药后病痛明显减轻的时间。
病人序号
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
用药剂量X1
2 2 2 2 2 2 5 5 5 5 5 5 7 7 7 7 7 7 10 10 10 10 10 10
性别X2 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1
血压组别X3
0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75
病痛减轻时间y
35 43 55 47 43 57 26 27 28 29 22 29 19 11 14 23 20 22 13 8 3 27 26 5
（要求对对所使用的命令加以说明） 解答：
基本模型：病痛减轻时间y与用药剂量x1，性别x2，血压组别x3之间的多元线性回归模型为
y??0??1x1??2x2??3x3??4x12?? （1）
其中?0,?1,?2,?3,?4是待估计的回归系数，?是随机误差。
模型求解：利用MATLAB统计工具箱中的命令regress求解：
x1=[2 2 2 2 2 2 5 5 5 5 5 5 7 7 7 7 7 7 10 10 10 10 10 10];x2=[0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1];
x3=[0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 ];
y=[35 43 55 47 43 57 26 27 28 29 22 29 19 11 14 23 20 22 13 8 3 27 26 5]';
22
x12=x1.^2;
x=[ones(24,1) x1' x2' x3' x12']; [b,bint,r,rint,stats]=regress(y,x); b,bint,stats 运行结果：
b = 63.1291 -10.2706 5.6667 -1.5000 0.5111 bint =48.7173 77.5409 -14.9243 -5.6169 -0.0213 11.3546 -15.4325 12.4325 0.1319 0.8903
stats = 0.8275 22.7903 0.0000 44.3109
结果分析：由结果可知，R?0.8275，即因变量的82.75%可由模型（1）确定，?2,?3的置信区间包含零点，因而模型从整体来看不理想，还有待改进。
进一步讨论：?2,?3的置信区间包含零点，说明基本模型（1）存在缺点，为寻找改进的方向，进行残差分析。将影响因素分为用药剂量与性别—血压组合两类，性别—血压组合的定义如表1。
表1 性别—血压组合
组合 性别 血压
1 0 0.25
2 1 0.25
3 0 0.5
4 1 0.5
5 0 0.75
6 1 0.75
2画出?与x1之间的关系，以及?与性别—血压之间的关系，对残差进行分析。
根据分析?与性别—血压之间的关系的点分布所以我们决定引入交互项，建立新的回归模型。 更好的模型：将对因变量y产生影响的所有交互项逐次加入模型，但只有x1x2变量的加入优化了模型。 将x1x3变量选入模型，建立如下回归模型
y??0??1x1??2x2??3x3??4x1x3??5x12 （2）
模型求解：利用MATLAB统计工具箱中的命令regress求解：
x1=[2 2 2 2 2 2 5 5 5 5 5 5 7 7 7 7 7 7 10 10 10 10 10 10]; x2=[0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1];
x3=[0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 ]; x1x3=x1.*x3; x12=x1.^2;
y=[35 43 55 47 43 57 26 27 28 29 22 29 19 11 14 23 20 22 13 8 3 27 26 5]'; x=[ones(24,1) x1' x2' x3' x1x3' x12']; [b,bint,r,rint,stats]=regress(y,x); b,bint,stats 运行结果：
b =40.5408 -6.5059 5.6667 43.6765 -7.5294 0.5111 bint = 26.8318 54.2499 -10.0338 -2.9780
23
1.8308 9.5025 22.1779 65.1751 -10.7522 -4.3066 0.2554 0.7668
stats = 0.9262 45.2100 0.0000 20.0014
模型分析：结果显示，R?0.9262指因变量（病痛减轻时间）的92.62%可由模型确定，p值为0，模型中的?0?40.5408，且它,?1??6.5059,?2?5.6667,?3?43.6765,?4??7.5294,?5?0.5111们的置信区间都不包含零点，因而模型从整体来看是可用的。
综上，模型为y?40.5408?6.5059x1?5.6667x2?43.6765x3?7.5294x1x3?0.5111x122。
预测结果：
x1=2:10,x2=0,x3=0.25: 得到病痛减轻时间：
y=36.7279 30.8953 26.0848 22.2966 19.5306 17.7868 17.0652 17.3658 18.6887
6．实验心得（质疑、建议）：
24