数学建模B作业全部-部分答案.doc

2015年数学建模B作业
（全部，共23题）

作业要求
1.作业解答写在实验报告纸上，无需抄题，但要写题号。
2.实验报告纸上要写程序，程序中可不抄数据。
3.将程序运行的重要结果有选择的展示在实验报告纸上，并做结果分析。
4.从第三周开始，每周要交1次作业。每次作业的题目根据进度由老师安排。如老师未作说明，那就是：课讲到哪里作业就做到哪里。
5.如何收作业，听任课老师安排。
6.不收作业的打印版、电子版。

第一部分 多元统计

2015-1 回归分析
某种水泥在凝固时放出的热量y(k/g)与水泥中的3CaOAl2O3的成分 （%），3CaOSiO2的成分x2（%）,4CaOAl2O3Fe2O3的成分x3（%）,2CaOSiO2的成分x4（%）的观测值如下表，试以y为因变量，以x1，x2，x3，x4为自变量建立多元回归方程并作显著性检验。

样本点
x1
x2
x3
x4
y
1
7
26
6
60
78.5
2
1
29
15
52
74.3
3
11
56
8
20
104.3
4
11
31
8
47
87.6
5
7
52
6
33
95.9
6
11
55
9
22
109.2
7
3
71
17
6
102.7
8
1
31
22
44
72.5
9
2
54
18
22
93.1
10
21
47
4
26
115.9
11
1
40
23
34
83.8
12
11
66
9
12
113.3
13
10
68
8
12
109.4
解：考虑到变量间可能存在多重共线性，采用逐步回归，程序如下：
data ex;input x1-x4 y @@;
cards;
7 26 6 60 78.5
1 29 15 52 74.3
11 56 8 20 104.3
11 31 8 47 87.6
7 52 6 33 95.9
11 55 9 22 109.2
3 71 17 6 102.7
1 31 22 44 72.5
2 54 18 22 93.1
21 47 4 26 115.9
1 40 23 34 83.8
11 66 9 12 113.3
10 68 8 12 109.4
;
proc reg;
model y=x1-x4/selection=stepwise;
run;
程序运行最终结果如下：

线性回归方程为：，由于几处pr>F的值均小于0.0001，故回归方程的线性性及各参数的显著性检验均通过，且拟合优度达到R2=0.9787。

2015-2 聚类分析
DNA是由A，T，C，G这4种碱基按一定顺序排成的序列，长短不一，其中碱基含量的百分比不同通常能揭示该序列的一些规律，试根据下表所给出的20条DNA序列的碱基含量百分比对其20条DNA序列进行分类。
（注，计算式下面的数据需要转置）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
A
33
30
30
47
26
39
39
31
23
20
39
36
28
33
32
40
39
32
24
22
T
15
17
7
32
12
14
21
21
17
15
55
55
57
55
71
51
29
55
62
62
C
19
18
24
12
26
14
11
18
23
30
5
3
11
9
0
9
27
13
16
19
G
44
46
50
20
47
44
40
41
48
45
11
16
14
13
7
10
15
10
8
7
解：依题意，首先将原始数据进行处理，转化成百分比，程序如下：
data ex;input a t c g@@;
cards;
0.2973 0.1351 0.1712 0.3964
0.2703 0.1532 0.1622 0.4144
0.2703 0.0631 0.2162 0.4505
0.4234 0.2883 0.1081 0.1802
0.2342 0.1081 0.2342 0.4234
0.3514 0.1261 0.1261 0.3964
0.3514 0.1892 0.0991 0.3604
0.2793 0.1892 0.1622 0.3694
0.2072 0.1532 0.2072 0.4324
0.1818 0.1364 0.2727 0.4091
0.3545 0.5000 0.0455 0.1000
0.3273 0.5000 0.0273 0.1455
0.2545 0.5182 0.1000 0.1273
0.3000 0.5000 0.0818 0.1182
0.2909 0.6455 0 0.0636
0.3636 0.4636 0.0818 0.0909
0.3545 0.2636 0.2455 0.1364
0.2909 0.5000 0.1182 0.0909
0.2182 0.5636 0.1455 0.0727
0.2000 0.5636 0.1727 0.0636
;
proc cluster method=single; /*最短距离法*/
proc tree;run;
聚类图如下，根据动态聚类图可以看出，此处20个DNA序列分成三类较为合适，具体情况如下：
第一类：4，17；
第二类：1，2，3，5，6，7，8，9，10；
第三类：11，12，13，14，15，16，18，19，20

2015-3 判别分析
观测3名健康人和4名心肌梗塞病人心电图的3项指标x1，x2，x3所得观测值如下表，试判别心电图3项指标为（400.72，49.46，2.25）的人属于两类中的哪一类，并指出哪个指标在判别分析中占有最重要的地位。
类
病人编号
x1
x2
x3
健康人
1
436.70
49.59
2.32
2
290.67
30.02
2.46
3
352.53
36.26
2.36
心肌梗塞病人
1
510.47
67.64
1.73
2
510.41
62.71
1.58
3
470.30
54.4.
1.68
4
364.12
46.26
2.09
解：根据判别分析编写程序如下：
data ex;input leibie x1 x2 x3@@;
cards;
1 436.70 49.59 2.32
1 290.67 30.02 2.46
1 352.53 36.26 2.36
2 510.47 67.64 1.73
2 510.41 62.71 1.58
2 470.30 54.4. 1.68
2 364.12 46.26 2.09
;
data ex1;input x1 x2 x3@@;
cards;
400.72 49.46 2.25
;
proc discrim data=ex testdata=ex1 anova manova simple list testout=ex2;
class leibie;
proc print data=ex2;run;
根据运行结果知，所建立的判别规则误判率为0：

最终判别结果为该人属于第一类：

2015-4 主成分分析
某市为全面分析机械类各企业的经济效益，选择了8个不同的利润指标，14个企业关于这8个指标的统计数据如下表，试进行主成分分析并将14个企业的经济效益进行排序。
企
业
净产值
利润率
固定资产
利润率
总产值
利润率
销售收入
利润率
产品成本
利润率
物耗利
润率
人均利
润率
流动资金
利润率
1
40.4
24.7
7.2
6.1
8.3
8.7
2.442
20.0
2
25.0
12.7
11.2
11.0
12.9
20.2
3.542
9.1
3
13.2
3.3
3.9
4.3
4.4
5.5
0.578
3.6
4
22.3
6.7
5.6
3.7
6.0
7.4
0.176
7.3
5
34.3
11.8
7.1
7.1
8.0
8.9
1.726
27.5
6
35.6
12.5
16.4
16.7
22.8
29.3
3.017
26.6
7
22.0
7.8
9.9
10.2
12.6
17.6
0.847
10.6
8
48.4
13.4
10.9
9.9
10.9
13.9
1.772
17.8
9
40.6
19.1
19.8
19.0
29.7
39.6
2.449
35.8
10
24.8
8.0
9.8
8.9
11.9
16.2
0.789
13.7
11
12.5
9.7
4.2
4.2
4.6
6.5
0.874
3.9
12
1.8
0.6
0.7
0.7
0.8
1.1
0.056
1.0
13
32.3
13.9
9.4
8.3
9.8
13.3
2.126
17.1
14
38.5
9.1
11.3
9.5
12.2
16.4
1.327
11.6

解：编写主成分分析的程序如下：
data ex;input x1-x8@@;
cards;
40.4 24.7 7.2 6.1 8.3 8.7 2.442 20.0
25.0 12.7 11.2 11.0 12.9 20.2 3.542 9.1
13.2 3.3 3.9 4.3 4.4 5.5 0.578 3.6
22.3 6.7 5.6 3.7 6.0 7.4 0.176 7.3
34.3 11.8 7.1 7.1 8.0 8.9 1.726 27.5
35.6 12.5 16.4 16.7 22.8 29.3 3.017 26.6
22.0 7.8 9.9 10.2 12.6 17.6 0.847 10.6
48.4 13.4 10.9 9.9 10.9 13.9 1.772 17.8
40.6 19.1 19.8 19.0 29.7 39.6 2.449 35.8
24.8 8.0 9.8 8.9 11.9 16.2 0.789 13.7
12.5 9.7 4.2 4.2 4.6 6.5 0.874 3.9
1.8 0.6 0.7 0.7 0.8 1.1 0.056 1.0
32.3 13.9 9.4 8.3 9.8 13.3 2.126 17.1
38.5 9.1 11.3 9.5 12.2 16.4 1.327 11.6
;
proc princomp out=prin;var x1-x8;
proc print data=prin;var Prin1-prin8;
run;
根据运行结果，以累积贡献率超过90%为标准，可选择三个主成分：

根据特征向量可以写出主成分表达式：

如第一主成分可写为如下，其它类似：

由变量前的系数大小可见，第一主成分主要是反映总产值利润率、销售收入利润率和产品成本利润率的，是用来衡量企业经营状况的一个综合指标，其它可类似分析。
另外，还可进行主成分得分分析，主成分得分的结果如下：

可见，在第一主成分上得分最高的是企业9，在第二主成分上得分最高的是企业1，在第三主成分上得分最高的是企业2。

2015-5 因子分析
有10例患者的4项肝功能指标的观测数据如下表，试作这4项指标的因子分析并对病人进行病情分析。

患者
转氨酶量
肝大指数
硫酸锌浊度
胎甲球
1
40
2.0
5
20
2
10
1.5
5
30
3
120
3.0
13
50
4
250
4.5
18
0
5
120
3.5
9
50
6
10
1.5
12
50
7
40
1.0
19
40
8
270
4.0
13
60
9
170
3.0
9
60
10
130
2.0
30
50
解：编写因子分析程序如下：
data ex;
input a b c d;
cards;
40 2.0 5 20
10 1.5 5 30
120 3.0 13 50
250 4.5 18 0
120 3.5 9 50
10 1.5 12 50
40 1.0 19 40
270 4.0 13 60
170 3.0 9 60
130 2.0 30 50
;
proc corr out=ex1;
proc factor data=ex1 outstat=ex2 method=prin priors=one rotate=orthomax score;
proc score data=ex score=ex2 out=ex3;
proc print;run;
根据程序结果，按累积贡献率超过90%，选择三个公因子：

为了便于解释，旋转过后的因子模式为：

由此可写出：，其它类似。
标准化因子得分系数如下：

由此有，其它类似。
根据上式有因子得分结果如下：

在三个公因子上得分最高的患者依次是：4，10，8。

2015-6 典型相关分析
棉花红铃虫第一代发蛾高峰日y1（元月1日到发蛾高峰日的天数）、第一代累计百株卵量y2、发蛾高峰日百株卵量y3及2月下旬到3月中旬的平均气温x1、1月下旬到3月上旬的日照小时累计数的常用对数x2的10组观测数据如下表，试作气象指标与虫情指标间的典型相关分析。

样本点
x1
x2
y1
y2
y3
1
9.2
2.01
186
46.3
14.3
2
9.1
2.2
169
30.7
14.0
3
8.6
2.3
171
144.6
69.3
4
10.2
2.2
171
69.2
22.7
5
5.6
2.1
181
16.0
7.3
6
6.1
2.2
174
2.7
1.3
7
8.2
2.1
172
26.3
7.9
8
8.8
1.9
186
247.1
85.2
9
9.7
2.1
176
53.6
25.3
10
10.3
2.2
161
62.7
29.3
解：编写程序如下：
data ex; input x1-x2 y1-y3@@;
cards;
9.2 2.01 186 46.3 14.3
9.1 2.2 169 30.7 14.0
8.6 2.3 171 144.6 69.3
10.2 2.2 171 69.2 22.7
5.6 2.1 181 16.0 7.3
6.1 2.2 174 2.7 1.3
8.2 2.1 172 26.3 7.9
8.8 1.9 186 247.1 85.2
9.7 2.1 176 53.6 25.3
10.3 2.2 161 62.7 29.3
;
proc cancorr;
var x1-x2;
with y1-y3;
run;
根据运行结果，按累计贡献率达到90%及显著性的要求，选定典型变量：

典型变量的表达式子可以根据以下结果写出：

，
可见V1主要反映的是x2,W1主要反映的是y3,由于V1和W1是关系最强的一对典型变量，也意味着x2对y3有着最重要的影响。从下面还可以看出各单个变量与典型变量间的关系强弱：

第二部分 非参数统计

2015-7 方法比较
某制造商想要比较两种不同的生产方法所花费的生产时间是否有差异。随机地选取了11个工人，每一个工人都分别使用两种不同的生产方法来完成一项相同的任务，在样本中的每一个工人都做了观察。数据见表，试用Wilcoxon秩和检验这两种方法有无差异？
工人编号
1
2
3
4
5
6
7
8
9
10
11
方法1
10.2
9.6
9.2
10.6
9.9
10.2
10.6
10.0
11.2
10.7
10.6
方法2
9.5
9.8
8.8
10.1
10.3
9.3
10.5
10.0
10.6
10.2
9.8
解：提出原假设，这两组方法没有显著性差异，用配对实验的符号检验法，相应代码如下：
data ex;
input x1 x2@@;
y=x1-x2;
cards;
10.2 9.5
9.6 9.8
9.2 8.8
10.6 10.1
9.9 10.3
10.2 9.3
10.6 10.5
10 10
11.2 10.6
10.7 10.2
10.6 9.8
;
proc univariate;
var y;
run;

运行结果如下：

从结果中可以看出，sign统计量为3，其显著性为0.1094，大于0.05，故接受原假设，认为这两组方法没有显著性差异。

2015-8 培训方案选择
为培训大学生志愿者为社区服务，设计了4种培训方案，记作为A,B,C,D.将报名的30名大学生随机地分为4组，分别接受不同培训。训练一周后，按规定的要求考试，评定的成绩如下，试用非参数检验方法检验这四种培训方案的有效性是否存在显著差异？
培训方案A
60,75,62,76,73,98,86
培训方案B
72,52,68,82,74,64,87
培训方案C
61,85,78,66,70,59,69,79
培训方案D
63,58,65,71,84,77,80,89
解：提出原假设，这四种培训方案方法没有显著性差异，相应代码如下：
data ex;
do a=1 to 4;input n@@;
do i=1 to n;
input x@@;
output;end;end;
cards;
7 60 75 62 76 73 98 86
7 72 52 68 82 74 64 87
8 61 85 78 66 70 59 69 79
8 63 58 65 71 84 77 80 89
;
proc npar1way wilcoxon;class a;var x;
run;
运行结果如下：

从结果中可以看出，Chi-Square统计量为0.5537，其显著性为0.9069，大于0.05，故接受原假设，认为四种培训方案方法没有显著性差异。

2015-9 双胞胎智力的相关分析
某研究所对10对双胞胎儿童的智力进行调查，试计算其Pearson、Spearman和Kendall相关系数并对其进行相关性检验。
双胞胎编号
先出生儿童X
后出生儿童Y
1
9.0
7.8
2
16.6
19.3
3
16.2
20.1
4
11.3
7.1
5
16.2
13.0
6
7.1
4.8
7
7.8
8.9
8
4.0
7.4
9
11.2
10.0
10
1.3
1.5

解：
求其Pearson，Spearman和Kendall相关系数，代码如下：
DATA new;
INPUT x y@@;
CARDS;
9.0 7.8
16.6 19.3
16.2 20.1
11.3 7.1
16.2 13.0
7.1 4.8
7.8 8.9
4.0 7.4
11.2 10.0
1.3 1.5
;
PROC CORR pearson spearman kendall;
VAR x y;
RUN;
结果如下：
Pearson Correlation Coefficients, N = 10
Prob > |r| under H0: Rho=0

x y

x 1.00000 0.88081
0.0008

y 0.88081 1.00000
0.0008

Spearman Correlation Coefficients, N = 10
Prob > |r| under H0: Rho=0

x y

x 1.00000 0.82067
0.0036

y 0.82067 1.00000
0.0036

Kendall Tau b Correlation Coefficients, N = 10
Prob > |r| under H0: Rho=0

x y

x 1.00000 0.67420
0.0071

y 0.67420 1.00000
0.0071
可见，x与y的Pearson相关系数为0.88081，概率为0.0008，达到极显著水平；Spearman相关系数为0.82067，概率为0.0036，达到极显著水平；Kendall相关系数0.67420，概率为0.0071达到极显著水平；故，x与y显著相关。

第三部分 预测预报

2015-10 灰色预测
陕西省农业总产值数据如下：
年份
1985
1986
1987
1888
1989
1990
1991
1992
1993
1994
总产值
62.9
58.8
61.4
87.2
104.9
124.8
110.7
129.0
155.3
219.03
请建立灰色系统GM（1，1）模型，并预测1995-1997三年的农业总产值。
解：有原始时间1985-1994序列，对生成1-AGO序列另外可得Yn见表：
、1-AGO序列、Yn
k
1
2
3
4
5
6
7
8
9
10

62.9
58.8
61.4
87.2
104.9
124.8
110.7
129
155.3
219.03

62.9
121.7
183.1
270.3
375.2
500
610.7
739.7
895
1114.03
Yn

58.8
61.4
87.2
104.9
124.8
110.7
129
155.3
219.03

利用MATLAB编程得：

function [X,c,error1,error2]=example9_11()
%利用MATLAB编程预测2003年中国蔬菜产量，
%并对预测结果做残差检验和后验差检验，程序如下：
X0=[62.9 58.8 61.4 87.2 104.9 124.8 110.7 129.0 155.3 219.03
];
k=3;
[X,c,error1,error2]=GM11(X0,k)
plot(1985:1994,X0,'g*-')
hold on
plot(1985:1997,X)
%%
function [X,c,error1,error2]=GM11(X0,k)
% 建立函数[X,c,error1,error2]=example9_3_2_3(X0,k)
% 其中X0为输入序列，k为预测长度，
% X为预测输出序列，c为后验差检验数，error1为残差，error2为相对误差
format long;
n=length(X0);
X1=[];
X1(1)=X0(1);
for i=2:n
X1(i)=X1(i-1)+X0(i); %计算累加生成序列
end
for i=1:n-1
B(i,1)=-0.5*(X1(i)+X1(i+1)); %计算B，Yn
B(i,2)=1;
Y(i)=X0(i+1);
end
alpha=(B'*B)^(-1)*B'*Y'; %做最小二乘估计
a=alpha(1,1);
b=alpha(2,1);
d=b/a; %计算时间响应函数参数
c=X1(1)-d;
X2(1)=X0(1);
X(1)=X0(1);
for i=1:n-1
X2(i+1)=c*exp(-a*i)+d;
X(i+1)=X2(i+1)-X2(i); %计算预测序列
end
for i=(n+1):(n+k)
X2(i)=c*exp(-a*(i-1))+d; %计算预测序列
X(i)=X2(i)-X2(i-1);
end
for i=1:n
error(i)=X(i)-X0(i);
error1(i)=abs(error(i)); %计算残差
error2(i)=error1(i)/X0(i); %计算相对误差
end
c=std(error1)/std(X0); %计算后验差检验数

运行结果见表格：

年份
1985
1986
1987
1888
1989
1990
1991
1992
1993
1994
1995
1996
1997
总产值
62.9
58.8
61.4
87.2
104.9
124.8
110.7
129
155.3
219.03

预测值
62.9
58.58326
68.15687
79.29499
92.25329
107.3292
124.8688
145.2748
169.0154
196.6357
228.7697
266.155
309.6498
残差
0
0.216
6.756
7.905
12.64
17.47
14.16
16.27
13.71
22.39

相对误差
0
0.003
0.110
0.090
0.120
0.139
0.127
0.126
0.088
0.102

画出预测与实际值变化曲线，如图所示：
预测与实际值变化曲线

实验模型以及结果检验：由表与图的结果可见，预测值与实际值偏离不大，其后验残差检验数C=0.1475小于0.35，所以模型精度为优。

2015-11 序列预测
某车站1993-1997年各月的列车运行数量数据如下表，试用时间序列建立合适的模型。并预测1998年1月的数值
1196.8 1181.3 1222.6 1229.3 1221.5 1148.4 1250.2 1174.4 1234.5 1209.7
1206.5 1204.0 1234.1 1146.0 1304.9 1221.9 1244.1 1194.4 1281.5 1277.3
1238.9 1267.5 1200.9 1245.5 1249.9 1220.1 1267.4 1182.3 1221.7 1178.1
1261.6 1274.5 1196.4 1222.6 1174.7 1212.6 1215.0 1191.0 1179.0 1224.0
1183.0 1288.0 1274.0 1218.0 1263.0 1205.0 1210.0 1243.0 1266.0 1200.0
1306.0 1209.0 1248.0 1208.0 1231.0 1244.0 1296.0 1221.0 1287.0 1191.0
解：
（1） 首先进行平稳性检验：
data a;/*a为数据名*/
input lieche@@;/*lieche为变量名*/
month=intnx('month','1jan1993'd,_n_-1);/*intnx间隔取时间变量*/
format month date.;/*月按？？？？*/
cards;
1196.8 1181.3 1222.6 1229.3 1221.5 1148.4 1250.2 1174.4 1234.5 1209.7
1206.5 1204.0 1234.1 1146.0 1304.9 1221.9 1244.1 1194.4 1281.5 1277.3
1238.9 1267.5 1200.9 1245.5 1249.9 1220.1 1267.4 1182.3 1221.7 1178.1
1261.6 1274.5 1196.4 1222.6 1174.7 1212.6 1215.0 1191.0 1179.0 1224.0
1183.0 1288.0 1274.0 1218.0 1263.0 1205.0 1210.0 1243.0 1266.0 1200.0
1306.0 1209.0 1248.0 1208.0 1231.0 1244.0 1296.0 1221.0 1287.0 1191.0
;
run ;
proc gplot;/*画图*/
plot lieche*month;/*纵轴为lieche,横轴为mouth*/
symbol v=square i=join c = red;/*图形特征，v表示点的形状，i表示图形连线的情况，c代表颜色*/
proc arima data = a;/*调用arima模块*/
identify var=lieche nlag = 22;/*延迟阶数为22阶*/
run;
运行得自相关图：

由此自相关图可看出，自相关系数很快的衰减向0，且始终控制在2倍范围内，可以认为该序列为平稳序列。
时序图：

由图可知，此车站列车运行数量数据在一个常数值附近随机波动，而且波动范围有界，无明显趋势及周期特征，基本可以视序列为平稳序列。

（2）进行随机性检验：
选取结果中The ARIMA Procedure部分：

由于统计量P值均大于0.05，则认为在0.05的显著水平下，无法拒绝原假设，即不能显著拒绝序列为纯随机序列的假定，因而认为此车站列车运行数量为纯随机波动序列，各序列之间没有任何行相关关系，即为无记忆序列，也就是说，该车站列车运行数量前后两年并无大的联系，也就是实说，我们很难根据历史信息预测未来年份此车站列车运行数量，故，该平稳序列不值得继续分析下去，对该序列分析到此结束。

2015-12 序列预测
对我国1952-1994年的社会消费品零售总额数据建立合适的时间序列模型，并预测1995-1997年的数据。
社会消费品零售总额
1952

262.7
328.8
356.1
1955
364.0
424.0
441.6
481.2
556.5
1960
595.4
537.7
543.7
544.8
572.7
1965
590.1
632.8
679.1
649.2
698.2
1970
728.8
776.9
853.5
917.7
967.4
1975
1046.4
1099.0
1174.3
1264.9
1476.0
1980
1794.0
2002.5
2181.5
2426.1
2899.2
1985
3801.4
4374.0
5115.0
6534.6
7074.2
1990
7250.3
8245.7
9704.8
12462.1
16264.7
解：（1）首先进行平稳性检验：
data a;/*a为数据名*/
input xf@@;/*xf为变量名*/
year=intnx('year','