止痛剂数学模型.doc

2013-2014-2数学建模实验综合实验

题 目 止痛剂回归统计模型

学 院
专业班级
姓 名
学 号

2014 年 6月 8日
止痛剂回归统计模型
一、 实验目的与意义：
1、练习实际问题的建模过程；熟悉数学建模步骤
2、练习Matlab基本编程命令；
二、 实验要求：
1、较能熟练应用Matlab基本命令和函数；
2、注重问题分析与模型建立，了解建模小论文的写作过程；
3、提高Matlab的编程应用技能。
三、 实验学时数：4学时
四、 实验类别：综合性
五、 实验内容与步骤：
问题 止痛剂是现在医学中常见的药剂，在世界范围内广泛使用的药剂，是人类健康存在必不可少的要素，因此有必要对止痛剂做一些研究。某一个医药公司的新药研究部门为了更好的掌握一种新止痛剂的疗效，设计了一个药物试验，给患同样病痛的病人使用这种新止痛剂的一下四个剂量中的某一个；2g,5g,7g和10g,并记录每个病人病痛明显减轻的时间（以分钟记），为了解新药的疗效与病人性别和血压有什么关系，试验过程中研究人员把病人按性别及血压的低、中、高三档平均分配进行测试。通过比较每个病人血压的历史数据，从高到低分成三组，分别记为0.25,0.50，0.75。试验结束后，公司的1结果见下表（性别以0表示女，1表示男），根据这些数据建立数学模型，来分析新药药效与病人性别和血压的关系，并进行统计预测。
表1 24名病人药剂量、性别、血压及疼痛时间数据
病人序号
病痛减轻时间/min
用药剂量/g
性别
血压组别
1
35
2
0
0.25
2
43
2
0
0.50
3
55
2
0
0.75
4
47
2
1
0.25
5
43
2
1
0.50
6
57
2
1
0.75
7
26
5
0
0.25
8
27
5
0
0.50
9
28
5
0
0.75
10
29
5
1
0.25
11
22
5
1
0.50
12
29
5
1
0.75
13
19
7
0
0.25
14
11
7
0
0.50
15
14
7
0
0.75
16
23
7
1
0.25
17
20
7
1
0.50
18
22
7
1
0.75
19
13
10
0
0.25
20
8
10
0
0.50
21
3
10
0
0.75
22
27
10
1
0.25
23
26
10
1
0.50
24
5
10
1
0.75
分析与假设 由于止痛剂在减轻病人的伤痛方面有很大帮助，病人关心疗效也是病痛减轻的时间长短，作为医药公司要对此有全面的了解，假设
（1）病痛减轻时间只与药剂量、性别、血压有关
（2）这24名病人是独立选取的，
（3）记x1(g)为药剂量，x2为性别，x3为血压组别，y（min)为病人病痛明显减轻的时间。
（4）为回归系数（i=1,2,3,4,5)
基本模型 为了大致地分析y与x1,x2,x3的关系，利用表1的数据分别做出y与x1,x2,x3的散点图（如下图）

图1 病痛减轻时间与药剂量的散点图

图2 病痛减少时间与性别的散点图

图3 病痛减少时间与血压组别的散点图
从图1可以发现，随着x1的增加，y有向下弯曲减少的趋势，故可考虑用二次函数模型

+ （1）
在图2,3中病痛减少时间与性别和血压无明显的关系，故不妨设其为线性关系，可考虑用线性模型
（2）
综合上面的分析，结合模型（1）和（2）建立如下的回归模型

（3）
（3）式右端的x1、x2和x3称为回归变量（自变量），是给定药剂量x1、性别x2和血压x3时，病痛减少时间y的平均值，其中，，，，称为回归系数，由表1的数据估计，影响y的其他因素作用都包含在随机误差中。如果模型选择的合适，应为服从均值为0的正态分布。
模型求解 直接利用matlab统计工具箱中命令regress求解，使用格式为
[b,bint,r,rint,stats]=regress(y,x,alpha)
其中输入y为模型（3）中y的数据（24维向量，n=24）,x为对应于回归系数的数据矩阵[1,x1,x2,x3,x1.^2](n*5矩阵，其中第一列为全1向量)，alpha为置信水平，b为的估计值，bint为b的置信区间，r为残差向量，rint为r的置信区间，stats为回归模拟的建议你统计量。
得到模型（3）的回归系数估计值及其置信区间（置信水平为0.05），结果见下表2
模型（3）的计算结果
参数
参数估计值
参数置信区间

63.1291
[48.7173 77.5409]

-10.2706
[-14.9243 -5.6169]

5.6667
[-0.0213 11.3546]

-1.5000
[-15.4325 12.4325]

0.5111
[0.1319 0.8903]
= 0.8275 F= 22.7903
P= 0.0000
= 44.3109
结果分析 表2显示，= 0.8275指因变量y的82.75%可由模型确定，p远远小于置信水平，因而模型（3）从整体来看是可用的。检查它们的置信区间发现，有，的置信区间包含零点，表明x2，x3不是太显著。
模型预测 将回归系数的估计值代入模型（3），
Y=63.1291-10.2706x1+5.6667x2-1.5000x3+0.5111x1^2 （4）
即可预测疼痛减少的时间, 设x1=5g，x2=0，x3=0.5，则疼痛减少时间y=63.1291-10.2706*5+5.6667*0-1.5000*0.5+0.5111*25=23.8036其置信度为95%的预测区间为[22.61342,24.99378]回归模型的一个重要应用是，对于给定的回归变量的取值，可以以一定的置信度预测因变量的取值范围，及预测区间。
模型改进 模型（3）中回归变量x1,x3对因变量y的影响是相互独立的，即疼痛减少时间y与药剂量的关系由回归系数确定，而不依赖于x3，同样y的均值与x3的线性关系由其回归系数确定，不依赖于x1，根据直觉和经验可以猜想，x2,x3的乘积代表它们的交互作用，于是将模型（3）增加一项，得到
+ （5）
下面我们用表1的数据估计模型（5）的系数，利用matlab的统计工具箱得到结果见下表3.
表3 模型（5）的计算结果
参数
参数估计
参数置信区间

40.5408
[26.8318 54.2490]

-6.5059
[-10.0338 -2.9780]

5.6667
[1.8308 9.5025]

43.6765
[22.1779 65.1751]

0.5111
[ 0.2554 0.7668]

-7.5294
[-10.7522 -4.3066]
= 0.9262 F= 45.2100
p= 0.0000
=20.0014
表3与表2的结果相比，有所提高，说明模型（5）比模型（3）有所改进。并且，所有参数的置信区间，特别是x1,x3的交互作用项的系数的置信区间不包含零点，所以有理由相信模型（5）比模型（3）更符合实际。
用模型（5）对疼痛减少时间做预测，设x1=5g，x2=0，x3=0.5，则疼痛减少时间y=40.5405-6.5059*5+5.6667*0+43.6765*0.5+0.5111*25-7.5294*2.5=23.8033其置信度为95%的预测区间为[22,613135,24.993465],与模型（3）结果相比，略有减少，预测区间较短
模型评价 共建立了两个模型，第二个比第一个模型要更符合实际，总的来说还是比较好的，可以作为预测模型，但在这两个模型中均未涉及x1和x2,x2和x3的交互项，虽然模型已经能很好的预测结果，但忽略了部分因素，故还有待完善的地方。模型有优点也有缺点，首先此模型为了解止痛剂的疗效提供了较好符合实际情况的数学模型，为疼痛时间的预测提供了简洁方便的工具，用简单的形式表示出事物之间复杂的关系，然而，由于自己数学知识和在计算机知识方面的匮乏，使得模型存在一定缺陷。
参考文献：
（1）《数学模型》（第四版）姜启源，谢金星，叶俊。高等教育出版社。
（2）《数学实验》（第二版）姜启源，张立平，何青，高丽。高等教育出版社。
（3）《概率论与数理统计》（第二版）程依明，濮晓龙，茆诗松
（4）《数学建模方法及应用》韩中庚，高等教育出版社。
（5）《数学模型建模分析》蔡常丰，科学出版社。
附录：
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.50 0.75 0.25 0.50 0.75 0.25 0.50 0.75 0.25 0.50 0.75 0.25 0.50 0.75 0.25 0.50 0.75 0.25 0.50 0.75 0.25 0.50 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];
plot(x1,y,'*')
plot(x2,y,'o')
plot(x3,y,'+')
x=1:1:24;
X=[ones(24,1) x1 x2 x3 x1.^2];
[b,brint,r,rint stats]=regress(y,X)
X=[ones(24,1) x1 x2 x3 x1.^2 x1.*x3];
[b,brint,r,rint stats]=regress(y,X)

运行结果b =

63.1291
-10.2706
5.6667
-1.5000
0.5111

brint =

48.7173 77.5409
-14.9243 -5.6169
-0.0213 11.3546
-15.4325 12.4325
0.1319 0.8903

r =

-9.2574
-0.8824
11.4926
-2.9240
-6.5490
7.8260
1.8211
3.1961
4.5711
-0.8456
-7.4706
-0.0956
3.0956
-4.5294
-1.1544
1.4289
-1.1961
1.1789
1.8407
-2.7843
-7.4093
10.1740
9.5490
-11.0760

rint =

-20.6395 2.1248
-13.6505 11.8858
0.6177 22.3676
-15.1031 9.2550
-18.9094 5.8113
-3.8146 19.4666
-11.0113 14.6534
-10.0627 16.4549
-8.0935 17.2356
-13.7026 12.0114
-20.3041 5.3629
-12.9593 12.7681
-9.6772 15.8684
-17.6934 8.6345
-14.0055 11.6967
-11.4155 14.2733
-14.5353 12.1431
-11.6717 14.0295
-10.3898 14.0712
-15.4859 9.9172
-19.1161 4.2975
-1.0161 21.3641
-2.3263 21.4244
-22.0557 -0.0963

stats =

0.8275 22.7903 0.0000 44.3109
b =

40.5408
-6.5059
5.6667
43.6765
0.5111
-7.5294

brint =

26.8318 54.2499
-10.0338 -2.9780
1.8308 9.5025
22.1779 65.1751
0.2554 0.7668
-10.7522 -4.3066

r =

-1.7279
-0.8824
3.9632
4.6054
-6.5490
0.2966
3.7034
3.1961
2.6887
1.0368
-7.4706
-1.9779
1.2132
-4.5294
0.7279
-0.4534
-1.1961
3.0613
-5.6887
-2.7843
0.1201
2.6446
9.5490
-3.5466

rint =

-9.2676 5.8117
-9.4997 7.7350
-3.3541 11.2805
-2.6136 11.8244
-14.5067 1.4087
-7.2928 7.8859
-4.7371 12.1440
-5.6742 12.0663
-5.8510 11.2285
-7.5960 9.6696
-15.6464 0.7052
-10.5680 6.6121
-7.4136 9.8401
-13.2476 4.1888
-7.9130 9.3689
-9.0993 8.1924
-10.1940 7.8018
-5.4458 11.5684
-12.7043 1.3269
-11.2959 5.7273
-7.4705 7.7107
-4.8257 10.1149
2.4226 16.6754
-10.9192 3.8260

stats =

0.9262 45.2100 0.0000 20.0014