【建模教程】-恒温箱温度变化的数学建模

----宋停云与您分享--------宋停云与您分享----全国第X届研究生数学建模竞赛题目 恒温箱温度变化的数学建模摘 要,对温度进行测量和控制也是科学实验和工业生产中经常需要解决的重要问题。随着计算机技术,,人们对高精度的温度测量和控制的要,其温度的测量和控制亦显得尤为重要和复杂。因此，本问题具有很强的研究背景和实际应用价值。为描述保温箱温度的变化规律，文章首先对环境温度进行三次样条插值，将环境温度的采样间隔变成和恒温箱一样，也就是一分钟。然后根据能量平衡原理和热力学6个未知量，取78156个方程式，用matlab求解这个超6个未知量的值，也就是恒温箱和隔热层的各项参数值。随后将剩余两个分钟的数据带入求得的方程组中进行检验，误差在接受范围内，说明建模合理。接下来，将热量平衡代数方程组变为功率平衡微分方程组，在考虑恒温箱控制条件的情况下用龙格库塔法进行求解，得到接下来恒温箱和隔热层的温度变换规律，画出温度变换曲线。最后对恒温箱的性能进行了评价。恒温箱隔热层热学平衡方程式超定方程组龙格库塔法----宋停云与您分享--------宋停云与您分享----参赛队号队员姓名陈* 管* 孙*（由组委会填写）----宋停云与您分享--------宋停云与您分享----XX大学承办----宋停云与您分享--------宋停云与您分享----恒温箱温度变化的数学建模一、问题重述,对温度进行测量和控制也是科学实验和工业生产中经常需要解决的重要问题。随着计算机技术,,人们对高精度的温度测量和控制的要,其温度的测量和控制亦显得尤为重要和复杂。因此，本问题具有很强的研究背景和实际应用价值。模型要求热源即加热器以恒定的速率加热恒温箱。恒温箱体内装有一个自动的温度控制器，在温度高于68.2华氏度会自动关闭加热系统,在温度低于67.8华氏度时会自动打开加热系统，但温度控制器读取温度有时间间隔，设定温度控制器读取温度有时间间隔为每1分钟1次。根据上述分析，本文拟解决如下问题：（1、根据已有的资料数据，利用所学知识及相关参考资料，通过研究温度传播规律,确定加热方案,（2、建立恒温箱体的温度变化模型。并根据小型试验数据及相关数学、物理方法，验证模型准确度和适用性；（3、对恒温箱产品的质量优劣做出评价；（4、对建立的模型进行误差分析，通过模型参数调节优化模型效果，并根据修改模型预测最有效数学模型。根据误差原因分析等对未来需要进行的试验和研究工作提出了一些建议；（5、根据相似准则设计相关实验，对已建立模型进行推广和扩充。以达到理论-实验-工程三者的互相验证。二、基本假设（1、环境温度不受热源、恒温箱和隔热层的影响；（2、忽略空气与恒温箱及隔热层之间的对流换热、热辐射；（3、所有的导热过程仅考虑稳态导热；（4、所有物体的物性不随温度变化；（567.8华氏度；（6、恒温箱体与箱顶隔热层温度始终不低于外界环境温度；----宋停云与您分享--------宋停云与您分享----三、符号约定c 恒温箱比热容1----宋停云与您分享--------宋停云与您分享----m1Tbox1

恒温箱质量t时间内恒温箱温度增量恒温箱箱体热导率----宋停云与您分享--------宋停云与您分享----Tbox

恒温箱某一时刻温度----宋停云与您分享--------宋停云与您分享--------宋停云与您分享--------宋停云与您分享----Tair

环境温度----宋停云与您分享--------宋停云与您分享----A 恒温箱与外界环境的接触表面积1 恒温箱体的厚度1----宋停云与您分享--------宋停云与您分享----Tgere

隔热层某时刻温度----宋停云与您分享--------宋停云与您分享----A 恒温箱体和隔热层的接触表面积2q 热源时间内所释放的热量c 隔热层的比热容2----宋停云与您分享--------宋停云与您分享----m2Tgere

隔热层的质量t时间隔热层温度的增量----宋停云与您分享--------宋停云与您分享----2 隔热层的热导率2 隔热层的厚度----宋停云与您分享--------宋停云与您分享----四、模型的建立与求解基本理论（1）热传导只要有温差,热量就会自发地从高温物体传到低温物体,所以传热是一种最常见的物,热是一种扩散效应,物体将能量通过热的方式释放出来,因此热其实是能量的传递。热的传递方式主要包括以下三种:热传导、热对流和热辐射。在对恒温箱的数学建模中，我们忽略热对流和热辐射，仅考虑热传导的主要作用。,在各部分之间没有发生相对宏观位移情况下,通过物质微观粒子()的热运动而进行的热量传;在气体和液体中,尽管也同样存在着导热现象。2恒温箱壁面导热示意图如图所示恒温箱壁，一块厚度为，横截面积为A的平板，两侧表面分别维持着均匀温度t 和t 实验表明单位时间内从表面1传递到表面2的热量Q与壁面间的温差w2tt t 以及导热面积A成正比，而与平板的厚度成反比，即w2Q=At方程中的比例系数反映了材料导热能力的大小，叫做导热系数或者热导率，单位是(mK)。导热率是一个物性参数，不同材料的导热率不同，即使是同一种材料，在不同温度下也具有不同的热导率。为了确保恒温箱内的光谱仪能够安全可靠的工作,我们需要将恒温箱内多余的热量3个主要方面。第一：组件级热设计，即对元器件级的热设计。第二：封装级设计，即对于电子控制模块、散热设备、PCB板的热设6----宋停云与您分享--------宋停云与您分享----（2）三次样条差值15min1min的数据。算法原理：----宋停云与您分享--------宋停云与您分享----设在区间[a,b]上给定n1个节点xi

(ax0

x x1 n

bxi

处的函数值----宋停云与您分享--------宋停云与您分享----为yf(x)(i0,1, ,n)。若函数S(x)满足如下三条：i i----宋停云与您分享--------宋停云与您分享----（1） 在每个子区间[xi1

,xi1,2, nS(x是三次多项式；i----宋停云与您分享--------宋停云与您分享----,n)；（2） S(x)y(,n)；i i（3） 在区间[abS(x的二阶导数Sx连续。S(xyf(x在区间[ab上的三次样条插值函数。----宋停云与您分享--------宋停云与您分享----子区间[xi1

,xi1,2, nS(x的表达式为：i----宋停云与您分享--------宋停云与您分享--------宋停云与您分享--------宋停云与您分享----S(x)

(xx)3i M

(xx )3i1 M

(y

h2 xxM)i iM)----宋停云与您分享--------宋停云与您分享----6h i1ih2

6hixx

i i1 6

i1 hi----宋停云与您分享--------宋停云与您分享----(yi

i M)6 i

i1, xh i1i

xxi----宋停云与您分享--------宋停云与您分享--------宋停云与您分享--------宋停云与您分享----关于参数Mi

的方程组（三弯矩方程组）：----宋停云与您分享--------宋停云与您分享----Mi1i

2Mi

Mi

i1

d, i1,2,,,n1----宋停云与您分享--------宋停云与您分享----h ih

, i1

1, hxx----宋停云与您分享--------宋停云与您分享----hi hhh

i hh

i i i1i----宋停云与您分享--------宋停云与您分享----i i1 i i16 y y yy----宋停云与您分享--------宋停云与您分享----d ( i1 i i i1)6f[x

,x,x ]----宋停云与您分享--------宋停云与您分享----i hhi i1

h hi1 i

i1

i1----宋停云与您分享--------宋停云与您分享----牛顿插值多项式：----宋停云与您分享--------宋停云与您分享----N(x)f[xn

]f[x0

,x](xx1 0

)f[x0

,x](xx1 0

)(xx)1----宋停云与您分享--------宋停云与您分享----f[x,x, ,x](xx)(xx

) (xx )----宋停云与您分享--------宋停云与您分享----0 1 k 0 1

k1----宋停云与您分享--------宋停云与您分享----f[x,x, ,x](xx)(xx0 1 n 0 1

) (xx ).n1----宋停云与您分享--------宋停云与您分享----构造牛顿插值多项式首先列出差商表，进而由差商表写出牛顿插值多项式。n阶差商为：----宋停云与您分享--------宋停云与您分享----f[x0

,x,,x,x]f[x,x,0 1,x]f[x,x,nnx x01,x ]n1n 0----宋停云与您分享--------宋停云与您分享----零阶差商和一阶差商：f[x]yf(x) f[x,x ]

f[x ]f[x]i1 i----宋停云与您分享--------宋停云与您分享----i i i

i i

x xi1 i----宋停云与您分享--------宋停云与您分享----（3）龙格库塔法算法原理：考虑用函数f(x,y)在若干点上的函数值的线性组合来构造近似公式，构造时要求近似公式在(xi,yi)TaylorxiTaylor展开式的前面几项重合，从而使近似公式达到所,在[xi,xi+1]这一步内多计算几个点的斜率值，然后将其进行加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格—库塔（Runge-Kutta）[xi,xi+1]上用四个点处的斜率加权平均作为平均斜率即局部截断误差是O(h5)。对于高阶常微分方程的初值问题：----宋停云与您分享--------宋停云与您分享----y(m)

f(x,y,y',y'', ,y(m1)), axb----宋停云与您分享--------宋停云与您分享----y(a)y,y'(a)y

', ,y(m1)(a)y

(m1),----宋停云与您分享--------宋停云与您分享----0 0 0将其转化为一阶微分方程组求解：y'y1 2y'y2 3y 'y,y)my,y)my', ,y(a)0 mm1 2y'f(x,y,m1 2----宋停云与您分享--------宋停云与您分享----y(a)y,y1 0 2

(a)

y(m1)0----宋停云与您分享--------宋停云与您分享----标准的四级四阶龙格—库塔法的向量形式：----宋停云与您分享--------宋停云与您分享----y y1(KK

K K)----宋停云与您分享--------宋停云与您分享----i1

i 6 1 2 3 4----宋停云与您分享--------宋停云与您分享----K hf(x,y)1 i i 1 1K hf(x h,y K) 2 i 2 i 2 1 1 1K hf(x2h,y2K)33 i i 2----宋停云与您分享--------宋停云与您分享----K4

hf(xi

h,yi

K)3----宋停云与您分享--------宋停云与您分享----其分量形式：----宋停云与您分享--------宋停云与您分享----y y 1(K K

K K )----宋停云与您分享--------宋停云与您分享---- j,i1

i,j 6 j1 j2

,y ,)nij,y ,)ni----宋停云与您分享--------宋停云与您分享----K hf(x;y,y ,j1 j i 2i,n h K K K,n----宋停云与您分享--------宋停云与您分享----K hf(x ;y

11,y

21, ,y

) j1,2,----宋停云与您分享--------宋停云与您分享---- j2 j i

2 2h K

2i 2K

ni 2K----宋停云与您分享--------宋停云与您分享----K hf(x2;y

12,y 22, ,y

n2)----宋停云与您分享--------宋停云与您分享---- j3 j i

2 2i 2

ni 2----宋停云与您分享--------宋停云与您分享----Kj4

hf(xijiy K113i113

,y K , ,y K )2i 23 ni n3----宋停云与您分享--------宋停云与您分享----（4）超定方程的最小二乘解小二乘法广泛地应用于工程计算中，用最小二乘法消除（平滑）误差，用最小二乘法从有噪声。甚至利用简单函数计算复杂函数的近似值，我们并不期望它的近似值多么精确（事实上很多时候也不用很精确，尽管如此还是希望计算出的近似数据与原始数据之间有相似之处。如果从线性代数角度来理解最小二乘法，实际上是将一个高维空间的向量投影到低维子空间所涉及的工作。超定方程组的最小二乘解ium×nGX=bGG=(g)m>n时即所谓的高矩阵，绝大多数情况下，超定方程组没r=b–GXium×n----宋停云与您分享--------宋停云与您分享----||bGX*||2超定方程组的最小二乘解是指正规方程组

min||bGX||XRm 2----宋停云与您分享--------宋停云与您分享----GTGX=GTb的解。如果系数矩阵(GTG)可逆，则正规方程组有唯一解。此时，最小二乘解可以形式地写为如下形式----宋停云与您分享--------宋停云与您分享----两种常用的方法如下

X=(GTG)-1GTb----宋停云与您分享--------宋停云与您分享----1．用对称矩阵的三角分解法解正规方程组GTGX=GTb；记A=GTG，则A是对称矩阵，由三角分解A=LDLT，其中L是下三角矩阵，D是对角矩阵。将这一算法写过三个过程：①解下三角方程组：LY1=GTb；Y1；=Y22QR分解直接求解超定方程组QR分解（正交三角分解QQR分解代入最小二乘解表达式中，得X=(RTQTQR)-1(QR)Tb=R-1QTb由此可知，用分解方法求超定方程组的最小二乘解只需求解上三角方程组RX=QTb----宋停云与您分享--------宋停云与您分享----模型168.2华氏度，恒温箱体温度高于箱顶隔热层，如图箱顶隔热层恒温箱体热源31恒温箱热量传播方向图模型268.267.8华氏度，热源停止加热；恒温箱体给箱顶隔热层加热，恒温箱及箱顶隔热层向外界空气散热，如图箱顶隔热层恒温箱体热源42恒温箱热量传播方向图模型3恒温箱及箱顶隔热层温度均高于67.8华氏度，热源停止加热；恒温箱温度低于箱顶隔热层，箱顶隔热层给恒温箱加热；恒温箱及箱顶隔热层向外界空气散热，如图----宋停云与您分享--------宋停云与您分享----箱顶隔热层箱顶隔热层恒温箱体热源53恒温箱热量传播方向图模型467.867.8华氏度，热源加热，热源与箱顶隔热层给恒温箱体加热，直到箱顶隔热层与恒温箱体温度相等。箱顶隔热层恒温箱体热源64恒温箱热量传播方向图模型建模15min测量一1min7所示。初始时68.2华氏度，恒温箱体温度高于箱顶隔热层。根据能量守恒，恒温箱体增加的热量与恒温箱对环境空气和箱顶隔热层的散热之和----宋停云与您分享--------宋停云与您分享----等于热源提供的热量。----宋停云与您分享--------宋停云与您分享----cmT

T )At/

T

)At/

qt （1）----宋停云与您分享--------宋停云与您分享----1 1

1 box

air 1

1 1 box

gere 2 1----宋停云与您分享--------宋停云与您分享----恒温箱体对箱顶隔热层的散热等于箱顶隔热层增加的热量与箱顶隔热层对环境空气的散热之和。----宋停云与您分享--------宋停云与您分享----cmT

T

)At/

(T

T )A

t/

（2）----宋停云与您分享--------宋停云与您分享----2 2 gere

1 box

gere 2

1 2 gere

air 2 2----宋停云与您分享--------宋停云与您分享----式中各变量含义分别为c1

是恒温箱比热容，m1

是恒温箱质量，Tbox

是时间内恒----宋停云与您分享--------宋停云与您分享----温箱温度增量，Tbox

是1

是恒温箱箱体热导率，T 是恒box----宋停云与您分享--------宋停云与您分享----温箱某一时刻温度，Tair

A是恒温箱与外界环境的接触表面积，1 1

是恒温----宋停云与您分享--------宋停云与您分享----箱体的厚度，Tgere

A2

是恒温箱体和隔热层的接触表面积，q是热----宋停云与您分享--------宋停云与您分享----源c2

m2

是隔热层的质量，Tgere

是时----宋停云与您分享--------宋停云与您分享----间隔热层温度的增量，2

是隔热层的厚度。2----宋停云与您分享--------宋停云与您分享----将上式进一步简化为AtBTbox

Cbox

Tair

)t(Tbox

Tgere

)t （3）----宋停云与您分享--------宋停云与您分享--------宋停云与您分享--------宋停云与您分享----DTgere

Egere

T

)t(Tbox

Tgere

)t （4）----宋停云与您分享--------宋停云与您分享--------宋停云与您分享--------宋停云与您分享----式中A

，B

cm1 11

，C

A，D

cm2 21

，E

1min----宋停云与您分享--------宋停云与您分享----121A A A A 1211 2 1 2 2 1 2 128078组数据，求解超定方程中的ABCDE五个参数。A=22.4495B=63.9245C=0.0366D=122.5048E=0.4669将能量守恒方程变换为功率守恒方程：----宋停云与您分享--------宋停云与您分享----TB boxtT

A(Tbox

Tgere

)C(Tbox

Tair

) （5）----宋停云与您分享--------宋停云与您分享----D t

(Tbox

Tgere

)E(Tgere

Tair

) ----宋停云与您分享--------宋停云与您分享--------宋停云与您分享--------宋停云与您分享----代入参数值可得：box63.9245dTboxdt

22.44951.0366Tbox

Tgere

0.0366Tair

（7）----宋停云与您分享--------宋停云与您分享----122.5048

dTgeredt

Tbox

1.4669Tgere

0.4669Tair

（8）----宋停云与您分享--------宋停云与您分享----79800.56240.31120.602，0.331，在接受范围内，验证了方程式的正确性。图7环境温度实际测量与三次采样根据式（7）和（8）80min的89所示。880min恒温箱温度拟合曲线----宋停云与您分享--------宋停云与您分享----980min箱顶隔热层温度拟合曲线23410。图10恒温箱与箱顶隔热层温度变化曲线----宋停云与您分享--------宋停云与您分享----图11恒温箱温度变化曲线68华氏度，箱顶隔热层与环境温度变化趋势相近。恒温箱体在刚开始持续受热，约在26068华氏度，然后在67.8华68.21011可知，随着环境温度的升高，恒恒温箱箱体温度趋近于稳定后（300分钟~1440分钟t317分钟，814分钟,317/814=0.4017。恒温箱箱体温度趋近于稳定后300分钟~1440分钟，最高温度68.4457华氏度，最低67.6212华氏度。模型的改进 恒温箱温度控制算法设计PID控制调节热源的加热速率。使恒PID控制算法简单、鲁棒性能好、可靠性高，被PID控制的原理框图。----宋停云与您分享--------宋停云与您分享----将给定值r与实际输出值c相减，得到控制偏差：etrcPID控制器输出的控制量，输入到被控对象进行控制：----宋停云与您分享--------宋停云与您分享----uK

e1tedtT

de----宋停云与您分享--------宋停云与您分享----其传递函数形式为：

p T 0I

D dt ----宋停云与您分享--------宋停云与您分享----1GsK 1 Ts1c p Ts pI----宋停云与您分享--------宋停云与您分享----式中，K

Tp I

是积分时间常数，TD

是微分时间常数。----宋停云与您分享--------宋停云与您分享----e（比例也不能过大，过大）积分环节：主要作用是消除静差，提高系统的控制精度。积分作用的强弱取决于积----宋停云与您分享--------宋停云与您分享----分时间常数TI

，其值越大，积分作用越弱；反之越强，引入积分调节可使系统稳定性----宋停云与您分享--------宋停云与您分享----下降，动态响应变慢。微分环节：能反应偏差信号的变化趋势（变化速率，并能在偏差信号变得太大之----宋停云与您分享--------宋停云与您分享----李嘉林.TEC制冷加热管制热的恒温箱设计[D].复旦大学,2012.陈国初,张琳,.恒温恒液位系统的动态机理建模与仿真[J].青岛大学学报(工程技术版)，2003,18(2):46-50.叶庆银,刘斌,臧润清等.小型恒温箱恒温性能的理论及实验研究 [J].制冷与空调,2007,7(4):48-50,81.DOI:10.3969/j.issn.1009-8402.2007.04.012.傅秦生.[M].北京：机械工业出版社，2007.李乃成，梅立泉.数值分析[M].北京：科学出版社，2011.附录三次样条插值代码：[hwx,grc]=textread('hwxwd.txt','%f%f','headerlines',1)wmwd=textread('wmwd.txt','%f')figureplot(t1,wmwd)T_interp=0:1:1440;wmwd1=interp1(t1,wmwd,T_interp,'linear');wmwd2=interp1(t1,wmwd,T_interp,'spline');plot(t1(1:7),wmwd(1:7),'o',T_interp(1:90),wmwd1(1:90),'r+',T_interp(1:90),wmwd2(1:90),'g*')温度曲线拟合代码：调用函数：Tair=wmwd2;wdqx(0,1000,34.6256,34.6256,Tair,1,1000)龙格库塔法数值积分函数：function[hwxwd,grcwd]=wdqx(x0,xn,Tbox0,Tgere0,Tair,h,n)%定义初值%h:步长；n:节点个数x=x0:h:xn;y(1,1)=Tbox0;%赋恒温箱、隔热层温度初值y(2,1)=Tgere0;i=1;whilei<=nwhile(i<=n)&&(y(1,i)<68.2)K1(:,1)=[h*func1(y(1,i),y(2,i),Tair(i)),h*func3(y(1,i),y(2,i),Tair(i))]';K2(:,1)=[h*func1(y(1,i)+K1(1,1)/2,y(2,i)+K1(2,1)/2,Tair(i)),h*func3(y(1,i)+K1(1,1)/2,y(2,i)+K1(2,1)/2,Ta----宋停云与您分享--------宋停云与您分享----ir(i))]';K3(:,1)=[h*func1(y(1,i)+K2(1,1)/2,y(2,i)+K2(2,1)/2,Tair(i)),h*func3(y(1,i)+K2(1,1)/2,y(2,i)+K2(2,1)/2,Tair(i))]';K4(:,1)=[h*func1(y(1,i)+K3(1,1),y(2,i)+K3(2,1),Tair(i)),h*func3(y(1,i)+K3(1,1),y(2,i)+K3(2,1),Tair(i))]';y(:,i+1)=y(:,i)+1/6*(K1+2*K2+2*K3+K4);%加热循环i=i+1;endwhile(i<=n)&&(y(1,i)>67.8)K1(:,1)=[h*func2(y(1,i),y(2,i),Tair(i)),h*func3(y(1,i),y(2,i),Tair(i))]';K2(:,1)=[h*func2(y(1,i)+K1(1,1)/2,y(2,i)+K1(2,1)/2,Tair(i)),h*func3(y(1,i)+K1(1,1)/2,y(2,i)+K1(2,1)/2,Tair(i))]';K3(:,1)=[h*func2(y(1,i)+K2(1,1)/2,y(2,i)+K2(2,1)/2,Tair(i)),h*func3(y(1,i)+K2(1,1)/2,y(2,i)+K2(2,1)/2,Tair(i))]';K4(:,1)=[h*func2(y(1,i)+K3(1,1),y(2,i)+K3(2,1),Tair(i)),h*func3(y(1,i)+K3(1,1),y(2,i)+K3(2,1),Tair(i))]';y(:,i+1)=y(:,i)+1/6*(K1+2*K2+2*K3+K4);%不加热循环i=i+1;endendhwxwd=y(1,:);grcwd=y(2,:);end微分方程1：functiondelta_hwx=func1(Tbox,Tgere,Tair)%加热方程时的恒温箱温度微分方程delta_hwx=1/63.9245*(22.4495-1.0366*Tbox+Tgere+0.0366*Tair);end微分方程2：functiondelta_hwx1=func2(Tbox,Tgere,Tair)%不加热方程时的恒温箱温度微分方程delta_hwx1=1/63.9245*(-1.0366*Tbox+Tgere+0.0366*Tair);end微分方程3：functiondelta_grc=func3(Tbox,Tgere,Tair)%隔热层温度微分方程----宋停云与您分享--------宋停云与您分享----delta_grc=1/122.5048*(Tbox-1.4669*Tgere+0.4669*Tair)end----宋停云与您分享--------宋停云与您分享----统计建模与R软件第五章习题答案(假设检验)Ex5.1>x<-c(220,188,162,230,145,160,238,188,247,113,126,245,164,231,256,183,190,158,224,175)>t.test(x,mu=225)OneSamplet-testdata: xt=-3.4783,df=19,p-value=0.002516alternativehypothesis:truemeanisnotequalto22595percentconfidenceinterval:172.3827211.9173sampleestimates:meanofx192.15原假设：油漆工人的血小板计数与正常成年男子无差异。备择假设：油漆工人的血小板计数与正常成年男子有差异。p值小于0.05，拒绝原假设，认为油漆工人的血小板计数与正常成年男子有差异。>t.test(x,mu=225,alternative="less")OneSamplet-testdata: xt=-3.4783,df=19,p-value=0.001258alternativehypothesis:truemeanislessthan22595percentconfidenceinterval:-Inf208.4806sampleestimates:meanofx192.15同样可得出油漆工人的血小板计数小于正常成年男子的结论。Ex5.2>pnorm(1000,mean(x),sd(x))[1]0.5087941>x[1]1067 9191196 7851126 936 9181156 920 948>pnorm(1000,mean(x),sd(x))----宋停云与您分享--------宋停云与您分享----[1]0.5087941x<=1000的概率为0.509,故x大于1000的概率为0.491.要点：pnorm计算正态分布的分布函数。在R软件中，计算值均为下分位点。Ex5.3>A<-c(113,120,138,120,100,118,138,123)>B<-c(138,116,125,136,110,132,130,110)>t.test(A,B,paired=TRUE)Pairedt-testdata: AandBt=-0.6513,df=7,p-value=0.5357alternativehypothesis:truedifferenceinmeansisnotequalto095percentconfidenceinterval:-15.62889 8.87889sampleestimates:meanofthedifferences-3.375p值大于0.05，接受原假设，两种方法治疗无差异。Ex5.4（1）正态性W检验：>x<-c(-0.7,-5.6,2,2.8,0.7,3.5,4,5.8,7.1,-0.5,2.5,-1.6,1.7,3,0.4,4.5,4.6,2.5,6,-1.4)>y<-c(3.7,6.5,5,5.2,0.8,0.2,0.6,3.4,6.6,-1.1,6,3.8,2,1.6,2,2.2,1.2,3.1,1.7,-2)>shapiro.test(x)Shapiro-Wilknormalitytestdata: xW=0.9699,p-value=0.7527>shapiro.test(y)Shapiro-Wilknormalitytestdata: yW=0.971,p-value=0.7754ks检验：>ks.test(x,"pnorm",mean(x),sd(x))----宋停云与您分享--------宋停云与您分享----One-sampleKolmogorov-Smirnovtestdata: xD=0.1065,p-value=0.977alternativehypothesis:two-sidedWarningmessage:Inks.test(x,"pnorm",mean(x),sd(x)):cannotcomputecorrectp-valueswithties>ks.test(y,"pnorm",mean(y),sd(y))One-sampleKolmogorov-Smirnovtestdata: yD=0.1197,p-value=0.9368alternativehypothesis:two-sidedWarningmessage:Inks.test(y,"pnorm",mean(y),sd(y)):cannotcomputecorrectp-valueswithtiespearson拟合优度检验，以x为例。>sort(x)[1]-5.6-1.6-1.4-0.7-0.5 0.4 0.7 1.7 2.0 2.5 2.5 2.8 3.0 3.5 4.0[16] 4.5 4.6 5.8 6.0 7.1>x1<-table(cut(x,br=c(-6,-3,0,3,6,9)))>p<-pnorm(c(-3,0,3,6,9),mean(x),sd(x))>p[1]0.048947120.249900090.620022880.900758560.98828138>p<-c(p[1],p[2]-p[1],p[3]-p[2],p[4]-p[3],1-p[4]);p[1]0.048947120.200952980.370122780.280735680.09924144>chisq.test(x1,p=p)Chi-squaredtestforgivenprobabilitiesdata: x1X-squared=0.5639,df=4,p-value=0.967Warningmessage:Inchisq.test(x1,p=p):Chi-squaredapproximationmaybeincorrectp值为0.967，接受原假设，x符合正态分布。（2）方差相同模型t检验：>t.test(x,y,var.equal=TRUE)----宋停云与您分享--------宋停云与您分享----TwoSamplet-testdata: xandyt=-0.6419,df=38,p-value=0.5248alternativehypothesis:truedifferenceinmeansisnotequalto095percentconfidenceinterval:-2.326179 1.206179sampleestimates:meanofxmeanofy2.065 2.625方差不同模型t检验：>t.test(x,y)WelchTwoSamplet-testdata: xandyt=-0.6419,df=36.086,p-value=0.525alternativehypothesis:truedifferenceinmeansisnotequalto095percentconfidenceinterval:-2.32926 1.20926sampleestimates:meanofxmeanofy2.065 2.625配对t检验：>t.test(x,y,paired=TRUE)Pairedt-testdata: xandyt=-0.6464,df=19,p-value=0.5257alternativehypothesis:truedifferenceinmeansisnotequalto095percentconfidenceinterval:-2.373146 1.253146sampleestimates:meanofthedifferences-0.56三种检验的结果都显示两组数据均值无差异。（3）方差检验：>var.test(x,y)Ftesttocomparetwovariances----宋停云与您分享--------宋停云与您分享----data: xandyF=1.5984,numdf=19,denomdf=19,p-value=0.3153alternativehypothesis:trueratioofvariancesisnotequalto195percentconfidenceinterval:0.63265054.0381795sampleestimates:ratioofvariances1.598361接受原假设，两组数据方差相同。Ex5.5>a<-c(126,125,136,128,123,138,142,116,110,108,115,140)>b<-c(162,172,177,170,175,152,157,159,160,162)正态性检验，采用ks检验：>ks.test(a,"pnorm",mean(a),sd(a))One-sampleKolmogorov-Smirnovtestdata: aD=0.1464,p-value=0.9266alternativehypothesis:two-sided>ks.test(b,"pnorm",mean(b),sd(b))One-sampleKolmogorov-Smirnovtestdata: bD=0.2222,p-value=0.707alternativehypothesis:two-sidedWarningmessage:Inks.test(b,"pnorm",mean(b),sd(b)):cannotcomputecorrectp-valueswithtiesa和b都服从正态分布。方差齐性检验：>Ftesttocomparetwovariancesdata: aandbF=1.9646,numdf=11,denomdf=9,p-value=0.3200alternativehypothesis:trueratioofvariancesisnotequalto195percentconfidenceinterval:----宋停云与您分享--------宋停云与您分享----0.50219437.0488630sampleestimates:ratioofvariances1.964622ab的方差相同。t检验：>t.test(a,b,var.equal=TRUE)TwoSamplet-testdata: aandbt=-8.8148,df=20,p-value=2.524e-08alternativehypothesis:truedifferenceinmeansisnotequalto095percentconfidenceinterval:-48.24975-29.78358sampleestimates:meanofxmeanofy125.5833 164.6000可认为两者有差别。Ex5.6二项分布总体的假设检验：>binom.test(57,400,p=0.147)Exactbinomialtestdata: 57and400numberofsuccesses=57,numberoftrials=400,p-value=0.8876alternativehypothesis:trueprobabilityofsuccessisnotequalto0.14795percentconfidenceinterval:0.10974770.1806511sampleestimates:probabilityofsuccess0.1425P值>0.05，故接受原假设，表示调查结果支持该市老年人口的看法Ex5.7二项分布总体的假设检验：>binom.test(178,328,p=0.5,alternative="greater")Exactbinomialtestdata: 178and328numberofsuccesses=178,numberoftrials=328,p-value=0.06794----宋停云与您分享--------宋停云与您分享----alternativehypothesis:trueprobabilityofsuccessisgreaterthan0.595percentconfidenceinterval:0.49576161.0000000sampleestimates:probabilityofsuccess0.5426829不能认为这种处理能增加母鸡的比例。Ex5.8利用pearson卡方检验是否符合特定分布：>chisq.test(c(315,101,108,32),p=c(9,3,3,1)/16)Chi-squaredtestforgivenprobabilitiesdata: c(315,101,108,32)X-squared=0.47,df=3,p-value=0.9254接受原假设，符合自由组合定律。Ex5.9利用pearson卡方检验是否符合泊松分布：>n<-length(z)>y<-c(92,68,28,11,1,0)>x<-0:5>q<-ppois(x,mean(rep(x,y)));n<-length(y)>p[1]<-q[1];p[n]=1-q[n-1]>chisq.test(y,p=p)Chi-squaredtestforgivenprobabilitiesdata: yX-squared=2.1596,df=5,p-value=0.8267Warningmessage:Inchisq.test(y,p=p):Chi-squaredapproximationmaybeincorrect重新分组，合并频数小于5的组：>z<-c(92,68,28,12)>n<-length(z);p<-p[1:n-1];p[n]<-1-q[n-1]>chisq.test(z,p=p)Chi-squaredtestforgivenprobabilitiesdata: zX-squared=0.9113,df=3,p-value=0.8227可认为数据服从泊松分布。----宋停云与您分享--------宋停云与您分享----Ex5.10ks检验两个分布是否相同：>x<-c(2.36,3.14,752,3.48,2.76,5.43,6.54,7.41)>y<-c(4.38,4.25,6.53,3.28,7.21,6.55)>ks.test(x,y)Two-sampleKolmogorov-Smirnovtestdata: xandyD=0.375,p-value=0.6374alternativehypothesis:two-sidedEx5.11列联数据的性检验：>x<-c(358,2492,229,2745)>dim(x)<-c(2,2)>chisq.test(x)Pearson'sChi-squaredtestwithYates'continuitycorrectiondata: xX-squared=37.4143,df=1,p-value=9.552e-10P值<0.05，拒绝原假设，有影响。Ex5.12列联数据的性检验：>y[,1][,2][,3][1,]451210[2,]462028[3,]282330[4,]111235>chisq.test(y)Pearson'sChi-squaredtestdata: yX-squared=40.401,df=6,p-value=3.799e-07P值<0.05，拒绝原假设，不，有关系。Ex5.13因有的格子的频数小于5，故采用fiser确切概率法检验性。>fisher.test(x)----宋停云与您分享--------宋停云与您分享----Fisher'sExactTestforCountDatadata: xp-value=0.6372alternativehypothesis:trueoddsratioisnotequalto195percentconfidenceinterval:0.046243825.13272210sampleestimates:oddsratio0.521271p值大于0.05，两变量，两种工艺对产品的质量没有影响。Ex5.14由于是在相同个体上的两次试验，故采用McNemar检验。>mcnemar.test(x)McNemar'sChi-squaredtestdata: xMcNemar'schi-squared=2.8561,df=3,p-value=0.4144p值大于0.05，不能认定两种方法测定结果不同。Ex5.15符号检验：H0：中位数>=14.6;H1:中位数<14.6>x<-c(13.32,13.06,14.02,11.86,13.58,13.77,13.51,14.42,14.44,15.43)>binom.test(sum(x)>14.6,length(x),al="l")Exactbinomialtestdata: sum(x)>14.6andlength(x)numberofsuccesses=1,numberoftrials=10,p-value=0.01074alternativehypothesis:trueprobabilityofsuccessislessthan0.595percentconfidenceinterval:0.00000000.3941633sampleestimates:probabilityofsuccess0.1拒绝原假设，中位数小于14.6Wilcoxon符号秩检验：>wilcox.test(x,mu=14.6,al="l",exact=F)----宋停云与您分享--------宋停云与您分享----Wilcoxonsignedranktestwithcontinuitycorrectiondata: xV=4.5,p-value=0.01087alternativehypothesis:truelocationislessthan14.6拒绝原假设，中位数小于14.6Ex5.16符号检验法：>x<-c(48,33,37.5,48,42.5,40,42,36,11.3,22,36,27.3,14.2,32.1,52,38,17.3,20,21,46.1)>y<-c(37,41,23.4,17,31.5,40,31,36,5.7,11.5,21,6.1,26.5,21.3,44.5,28,22.6,20,11,22.3)>binom.test(sum(x>y),length(x))Exactbinomialtestdata: sum(x>y)andlength(x)numberofsuccesses=14,numberoftrials=20,p-value=0.1153alternativehypothesis:trueprobabilityofsuccessisnotequalto0.595percentconfidenceinterval:0.45721080.8810684sampleestimates:probabilityofsuccess0.7接受原