第六讲概率统计建模

概率统计建模方法及其在Matlab中的实现

第六讲2002B彩票中的问题（古典概型）2004B电力市场的输电阻塞管理（多元线性回归）2005A长江水质的评价和预测（多元统计综合评价）2005BDVD在线租赁（概率分布-正态分布等）2006A出版设资源管理配置（主成分分析、方差分析）2006B艾滋病的评价体系及疗效的预测（统计回归分析）历年的建模竞赛题统计概率分布（古典概型、二项分布、正态分布等）随机服务模型（排队服务模型）时间序列模型（马氏模型）回归模型（一元、多元、逐步回归）聚类分析（主成分分析、方差分析）常用的概率模型SPSS、Matlab一、样本总体1、总体：人们研究对象的全体。2、个体：总体中的每一个基本单位。3、样本：从总体中随机产生的若干个个体的集合。统计的主要任务从样本推断总体二、频数表和直方图1、频数：将数据的取值范围划分为若干个区间，统计这组数据在每个区间出现的次数。2、直方图：以数据的取值为横坐标，频数为纵坐标画出的阶梯形图。区间的划分有等距划分和非等距划分。3、直方图的matlab实现（1）数据输入：直接输入——针对数据较少间接输入——针对数据较多（先写一个纯文本数据）例学生的身高和体重学校随机抽取100名学生，测量他们的身高和体重，所得数据如下表身高体重身高体重身高体重身高体重身高体重1727517162……17764

5516867……184701696416552……166491716516962……171711674716865……1705920行100名学生的身高和体重表先把上面表格里的数据保存在txt文本里，再在matlab里面导入该数据；最后整理数据。high=data(:,1:2:9);high=high(:);weight=data(:,2:2:10);weight=weight(:);（2）直方图命令：[N,X]=hist(Y,M)样本数组将[min(Y),max(Y)]等分为M份，缺省时默认为10返回M个小区间的中点返回M个小区间的频数Matlab程序loaddata.txt;high=data(:,1:2:9);high=high(:);weight=data(:,2:2:10);weight=weight(:);[n1,x1]=hist(high)[n2,x2]=hist(weight)hist(high)hist(weight)三、统计量统计量：反映样本数量特征的函数，它不含任何未知量。1、算术平均值和中位数——表示位置的统计量平均值：命令形式：mean(x)功能：返回x的均值命令形式：nanmean(x)功能：返回除了NaN外x的均值中位数：将数据从小到大排列后位于中间位置的数。命令形式：median(x)功能：返回x的中位数命令形式：nanmedian(x)功能：返回除了NaN外x的中位数三、统计量统计量：反映样本数量特征的函数，它不含任何未知量。2、标准差、方差和极差——表示变异程度的统计量标准差：命令形式：std(x)功能：返回x的标准差命令形式：std(x,1)功能：返回命令形式：var(x)功能：返回x的方差命令形式：var(x,1)功能：返回std(x,1)的平方方差：标准差的平方。极差：最大值与最小值之差。命令形式：range(x)功能：返回x的极差三、统计量统计量：反映样本数量特征的函数，它不含任何未知量。3、中心矩、偏度和峰度——表示分布形状的统计量中心矩：命令形式：moment(x,n)功能：返回x的n阶矩。注：偏度反映分布的对称性。V>0时，称为右偏度，此时数据位于均值右边比位于左边的多。V<0，称为左偏度。V接近于0，则认为分布时对称的。如正态分布，V=0。偏度：随机变量的标准化的3阶中心矩命令形式：skewness(x)功能：返回x的偏度。注：峰度时分布形状的另一种度量。正态分布的峰度是3。若V比3大的多，表示分布有沉重的尾巴，说明样本中有较多远离均值的数据。因而峰度可以作为衡量偏离正态分布的尺度之一。峰度：随机变量的标准化的4阶中心矩命令形式：kurtosis(x)功能：返回x的峰度。四、常见的概率分布1、正态分布：命令形式1：normpdf(x,u,v)功能：计算参数为u,v的正态分布密度函数在x处的值命令形式2：normcdf(x,u,v)功能：计算参数为u,v的正态分布的累积分布函数的值命令形式3：norminv(a,u,v)功能：计算临界值x命令形式4：normrnd(u,v,m,n)功能：产生服从参数为u,v的正态分布的mn的矩阵注：再matlab工具箱中，对每一种分布都提供了如下的几类函数。pdf—概率密度cdf—分布函数inv—逆累积分布函数rnd—随机数生成注：求某个分布的概率密度、分布函数等的格式：

namepdf(),namecdf()……或者pdf(‘name’,…),cdf(‘name’,…)…….常见的分布函数表name函数说明name函数说明betaBeta分布bino二项分布exp指数分布geo几何分布hyge超几何分布poiss泊松分布unif均匀分布unid离散均匀分布chiz卡方分布fＦ分布gamGamma分布norm正态分布tT分布logn对数正态分布nbin负二项分布ncf非中心Ｆ分布nct非中心t分布ncx2非中心卡方分布rayl瑞利分布weibWeibull分布五、参数估计1、点估计：用样本统计量确定总体参数的值。它是用一个值去估计另一个值，所以称为点估计。2、区间估计：称为的置信区间为置信水平为显著水平命令形式：[]=namefit(x,alpha)样本数据显著水平，缺省时默认为0.05例如[mu,sigma,muci,sigmaci]=normfit(x,alpha)返回均值u的点估计返回标准差v的点估计这两者的区间估计例：分别用金球、铂球测定引力常数(2)用金球测定观测值为：6.683,6.681,6.676,6.678,6.679,6.672;(2)用铂球测定观测值为：6.661,6.661,6.667,6.667,6.664;设测定总体服从正态分布，其参数未知，分别求该参数的置信度为0.9的置信区间。X=[6.683,6.681,6.676,6.678,6.679,6.672];Y=[6.661,6.661,6.667,6.667,6.664];[a1,b1,c1,d1]=normfit(X,0.1)[a2,b2,c2,d2]=normfit(Y,0.1)六、假设检验1、单个总体均值u的检验

原假设为：备选假设：（1）已知，关于u的检验命令形式：

[h,p,c]=ztest(x,mu,sigma,alpha,taic)接受与否的参数在原假设条件下样本均值出现的概率均值的置信区间样本均值标准差显著水平备选假设的选择注：（2）未知，关于u的检验命令形式：

[h,p,c]=ttest(x,mu,alpha,taic)例某种电子元件的寿命x（以小时计）服从正态分布，其均值和方差均未知。现得16只元件的寿命如下：280101212224379179264222362168250149260485170问：是否有理由认为该元件的平均寿命大于225小时？x=[159280101212224379179264222362168250149260485170];[h,p,c]=ttest(x,225,0.05,1)结论：拒绝原假设，认为寿命不大于225小时。2、双正态总体均值的假设检验比较两个方差相等的正态总体的均值是否相等(T检验)命令格式：[H,P,ci,stats]=ttest2(X,Y,alpha,tail)功能：对两个正态分布总体的采样X、Y进行T检验，对H,P,alpha的解释同上；tail是假设的备选项(即备择假设)，有三个值：tail=0是默认值，可省略,说明备选项为“均值不相等”；tail=1,说明备选项为“X的均值大于Y的均值”；tail=-1,说明备选项为“X的均值小于Y的均值”。ci给出均值差的置信区间；stats是个结构，包含以下元素：tstat(统计值)、df(自由度)。例

某灯泡厂在采用一项新工艺前后，分别抽取了10只进行寿命试验，寿命分别为：旧灯泡：2461,2404,2407,2439,2394,2401,2543,2463,2392,2458新灯泡：2496,2485,2538,2596,2556,2582,2494,2528,2537,2492假设灯泡的寿命服从正态分布，能否认为采用新工艺后，灯泡的寿命提高了？(a=0.01)x=[2461,2404,2407,2439,2394,2401,2543,2463,2392,2458];y=[2496,2485,2538,2596,2556,2582,2494,2528,2537,2492];alpha=0.01;[h,p,ci,st]=ttest2(x,y,alpha,-1)结果：h=1%拒绝原假设即认为寿命未提高p=6.3361e-005%p很小，对假设置疑ci=-Inf-44.6944st=tstat:-4.8567df:183、两个总体一致性的假设检验比较两个不知道确切分布的总体均值是否相等命令格式：[P,H,stats]=ranksum(X,Y,alpha)功能：对两个总体的采样X、Y进行检验，对H,P,alpha的解释同上；stats是个结构，包含二个元素：zval(均值差的正态统计值)和ranksum(统计的秩和值)。例

两台机床加工同一种轴，抽样测量产品的直径(mm):机床甲：33.592,33.862,33.751,33.673,33.847,33.778,33.631,33.911,33.785,33.928机床乙：34.221,33.947,33.856,34.039,34.000,33.924,34.125,34.273,33.968,33.923在a=0.05下能否认为两台机床加工的直径没有显著不同？clear;x=[33.592,33.862,33.751,33.673,33.847,33.778,33.631,33.911,33.785,33.928];y=[34.221,33.947,33.856,34.039,34.000,33.924,34.125,34.273,33.968,33.923];alpha=0.05;[p,h,st]=ranksum(x,y,alpha)结果：p=7.6854e-004%p很小，对假设置疑h=1%拒绝原假设即认为直径没有显著不同st=zval:-3.3639ranksum:604、两个样本具有相同连续分布的假设检验检验两个样本是否具有相同的连续分布命令格式：

[H,P,ksstat]=kstest2(X,Y,alpha,tail)功能：对两个总体的采样X、Y进行检验，对H,P,alpha的解释同上；tail是假设的备选项(即备择假设)，有三个值：tail=0是默认值，可省略,说明备选项为“不相等”；tail=1,说明备选项为“大于”；tail=-1,说明备选项为“小于”。ksstat表示测试统计量的值。clear;x=randn(1,10);y=randn(1,10)+x;[h,p,kst]=kstest2(x,y)例

两个正态分布的检验结果：h=0%接受原假设即认为两样本具有相同类型的连续分布p=0.6751

%表示假设成立的概率为0.6751kst=0.3000

5、正态分布的假设检验检验样本是否具有某种连续分布命令格式1：

[H,P,jbstat,cv]=jbtest(X,alpha)功能：对采样X进行检验是否服从正态分布，对H,P,alpha的解释同上；jbstat表示测试统计量的值;cv为是否拒绝假设的临界值。适合大样本命令格式2：

[H,P,lstat,cv]=lillietest(X,alpha)功能：对采样X进行检验是否服从正态分布，对H,P,alpha的解释同上；jbstat表示测试统计量的值;cv为是否拒绝假设的临界值。适合小样本。clear;m1=ones(1,11)*2.55;m2=ones(1,12)*2.65;m3=ones(1,17)*2.75;m4=ones(1,19)*2.85;m5=ones(1,26)*2.95;m6=ones(1,24)*3.05;m7=ones(1,22)*3.15;m8=ones(1,19)*3.25;m9=ones(1,13)*3.35;M=[m1,m2,m3,m4,m5,m6,m7,m8,m9];[h,p,lst,cv]=lillietest(M)hist(M)例从一批零件中随机抽取一组样品，下面是零件样品直径的统计表。在显著水平a=0.05下能否认为这批零件的直径服从正态分布？绘出统计数据的直方图。直径2.552.652.752.852.953.053.153.253.35频数111217192624221913结果：h=1%拒绝原假设即认为直径不服从正态分布p=Nan

%表示假设成立的概率很小lst=0.1062cv=0.0694%测试统计值大于临界值也表明应拒绝hist(M,n)---绘制向量M的直方图,n定义条方的数目,默认为106、概率分布的假设检验命令格式：

[H,P,ksstat,cv]=kstest(X,cdf,alpha,tail)功能：对采样X进行检验是否服从名为cdf类型的连续累积概率分布，cdf缺省为[],默认为标准正态分布,声明格式为两个相同长度的列向量：采样和采样对应的分布函数；对H,P,alpha，ksstat,cv的解释同上。例

clear;mu=1;sigma=2;x=normrnd(mu,sigma,20,1);alpha=0.01;lbd=3;[h,p,ksst,cv]=kstest(x,[x,expcdf(x,lbd)],alpha,0)%检验是否符合参数为3的指数分布q-q图：用qqplot函数生成两个样本的q-q（quan-tile分位数）图。若两样本来自同一分布，图中数据点呈直线关系，否则为曲线关系。qqplot(X,Y)：显示X和Y两个样本的q-q图。qqplot(X)：显示X的样本值与服从正态分布的理论数据之间的q-q图。例x=normrnd(0,1,100,1);y=normrnd(0.5,2,50,1);z=weibrnd(2,0.5,100,1);subplot(2,2,1),qqplot(x)holdon,subplot(2,2,2)qqplot(x,y),holdonsubplot(2,2,3),qqplot(z)holdon,subplot(2,2,4)qqplot(x,z)holdoff由第一个子图看出X服从正态分布。由第二个子图看出X和Y可看作同分布的。由第三个子图看出Z不服从正态分布。由第四个子图看出X和Z不是同分布的。1、单因素方差分析模型：

六、方差分析用来对比因变量在不同组中的平均值的统计方法单因素方差分析命令形式命令格式：[P,anovatab,stats]=anova1(X,group,displayopt)功能：比较多组数据的均值，返回这些均值相等的概率，从而判断因素对结果是否有显著影响。X为输入数据，列向量表示相互独立的样本观测值，具有相同长度；P为X的各列均值相等的概率，P越小，则质疑原假设（即均值不相等），表示因素的影响显著；group是与X对应的字符或字符串数组，用来声明X每一列中数据的名字或意义，可以省略；displayopt表示参数：on表示显示图，off表示隐藏图；anovatab返回方差分析表；stats返回一个附加的统计数据结构。例

将同一批同种牌号丝袜在不同温度下作弹力试验，得到数据表：温度试验30º40º50º60º70º80º14.36.110.06.59.39.527.87.34.88.38.78.833.24.25.48.67.211.446.54.19.68.210.17.8试检验温度对弹力有无显著影响。(α=0.05)clear;X=[4.3,6.1,10.0,6.5,9.3,9.5;7.8,7.3,4.8,8.3,8.7,8.8;3.2,4.2,5.4,8.6,7.2,11.4;6.5,4.1,9.6,8.2,10.1,7.8];[p,tab,stats]=anova1(X,[],'on')结果：p=0.0214%p很小，拒绝原假设tab='Source''SS''df''MS''F''Prob>F''Columns'55.5471511.10943.5254[0.0214]'Error'56.7225183.1512[][]'Total'112.269623[][][]stats=gnames:[6x1char]n:[444444]source:'anova1'means:[5.45005.42507.45007.90008.82509.3750]

df:18s:1.7752%总体标准差的无偏估计例将四种工艺下生产的灯泡进行寿命测试，得到数据表：工艺试验A1A2A3A4116201580146015002167016001540155031700164016201610417501720168051800试检验工艺对寿命有无显著影响。(α=0.05)clear;X=[1620,1670,1700,1750,1800,1580,1600,1640,1720,1460,1540,1620,1500,1550,1610,1680];group=[1,1,1,1,1,2,2,2,2,3,3,3,4,4,4,4];[p,tab,stats]=anova1(X,group,'on')结果：p=0.0331tab='Source''SS''df''MS''F''Prob>F''Groups'[62820][3][20940][4.0608][0.0331]'Error'[61880][12][5.1567e+003][][]'Total'[124700][15][][][]stats=gnames:{4x1cell}n:[5434]source:'anova1'means:[1708163515401585]df:12s:71.80992、双因素方差分析模型：双因素方差分析命令格式命令格式：

[P,anovatab]=anova2(X,reps,displayopt)功能：判断因素对结果是否有显著影响。X为输入数据，列向量表示因素1的差异，行向量表示2的差异；P是概率向量，P越小，则质疑原假设，表示因素的影响显著；reps声明每一状态下的试验次数；displayopt表示参数：on表示显示图，off表示隐藏图；anovatab返回方差分析表。例有3个工人分别在4台机器上加工某种零件，工作的3天中日产量列表如下：

工人B机器AB1B2B3A1151517191916161821A21717

171515

15192222A31517161817161818

18A41820221516171717

17试检验操作工人的技术水平有无显著差异？机器性能有无显著差异？交互作用的影响是否显著？(α=0.05)clear;A1=[15,15,17,19,19,16,16,18,21];A2=[17,17,17,15,15,15,19,22,22];A3=[15,17,16,18,17,16,18,18,18];A4=[18,20,22,15,16,17,17,17,17];X=[A1',A2',A3',A4'];reps=3;[p,tab]=anova2(X,reps,'on')解释：操作工人的技术水平无显著差异p=0.6645机器性能有显著差异p=0.0023交互作用的影响显著p=0.0002直线拟合：a=polyfit(x,y,1),b=polyfit(x,z,1),同一条直线y=0.33x+0.96(z=0.33x+0.96)从拟合到回归x=[01234],y=[1.01.31.52.02.3](+号)x=[01234],z=[0.61.950.92.851.8]（*号）问题：你相信哪个拟合结果？怎样给以定量评价?得到a=0.330.96b=0.330.96七、回归分析一元线性回归分析的主要任务是：模型参数估计1、回归系数的最小二乘估计

其中åå====niiniiynyxnx111,1，åå====niiiniiyxnxyxnx11221,1.

检验、预测与控制1、回归方程的显著性检验（Ⅰ）F检验法

（Ⅱ）t检验法（Ⅲ）r检验法2、回归系数的置信区间3、预测与控制（1）预测（2）控制收集一组包含因变量和自变量的数据；选定因变量与自变量之间的模型，利用数据按照最小二乘准则计算模型中的系数；利用统计分析方法对不同的模型进行比较，找出与数据拟合得最好的模型；判断得到的模型是否适合于这组数据,诊断有无不适合回归模型的异常数据；利用模型对因变量作出预测或解释。回归分析的主要步骤

为了研究钢材消费量与国民收入之间的关系，在统计年鉴上查得一组历史数据。例：钢材消费量与国民收入的关系

年份196419651966……197819791980消费(吨)698872988……144627362825收入(亿)109712841502……294831553372

试分析预测若1981年到1985年我国国民收入以4.5%的速度递增，钢材消费量将达到什么样的水平？

钢材消费量--------试验指标(因变量)Y；国民收入-----------自变量x；建立数据拟合函数y=E（Y|x）=f(x)；作拟合曲线图形分析。

问题分析：多元线性回归

b=regress(Y,X)1、确定回归系数的点估计值：MATLAB统计工具箱常用命令3、画出残差及其置信区间：

rcoplot（r，rint）2、求回归系数的点估计和区间估计、并检验回归模型：

[b,bint,r,rint,stats]=regress(Y,X,alpha)回归系数的区间估计残差用于检验回归模型的统计量，有三个数值：相关系数r2、F值、与F对应的概率p置信区间显著性水平（缺省时为0.05）rcoplot(r,rint)残差及其置信区间作图MATLAB7.0版本s增加一个统计量:剩余方差s2.使用命令regress实现一元线性回归模型的计算

b=regress(Y,X)或

[b,bint,r,rint,stats]=regress(Y,X,alpha)回归系数beta以及它们的置信区间残差向量r=Y-Y及它们的置信区间相关系数R2，F-统计量和与F(1,n-2)分布大于F值的概率p，p<时回归模型有效.默认值是0.05模型求解输入：(hg1.m)x=[10971284150213941303155519172051211122862311200324352625294831553372];y=[698872988807738102513161539156117651762196019022013244627362825];X=[ones(size(x')),x'],pause[c,cint,r,rint,stats]=regress(y',X,0.05),pausercoplot(r,rint)输出：c=-460.5282(参数a)0.9840(参数b)cint=-691.8478-229.2085(a的置信区间)0.87791.0900(b的置信区间)r=[79.124869.1244-29.3788-104.1112-83.5709-44.5286-109.7219-18.5724-55.6100-23.8029-51.4019449.6576-33.4128-109.36515.816092.1364-32.3827]’(残差向量)rint=（略）（参见残差分析图）stats=0.9631(R2)391.2713(F)0.0000(P{χ0})第12个数据点异常，可删出预测x1(1)=3372;(hgy1.m)fori=1:5x1(i+1)=1.045*x1(i);%未来五年国民收入以4.5%的速度递增

y1(i+1)=-460.5282+0.9840*x1(i+1);%钢材的预测值endx1,y1结果x1=3372.03523.73682.33848.04021.24202.1y1=3006.83162.93325.93496.33674.4变量选择影响因变量的因素：自变量x1,x2,xm及其简单函数,如

将所有影响显著的因素都纳入回归模型；最终的模型尽量简单,即包含尽量少的因素。变量选择的标准

从候选集合S={x1,…xk}中选出一子集S1(含pk个自变量)与因变量y构造回归模型,其优劣由s2度量.影响显著的自变量进入模型时，Q明显下降，s减小；影响很小的自变量进入模型时，Q下降不大，p的增加会使s变大.变量选择与逐步回归

逐步回归从候选集合中确定一初始子集；从子集外（候选集合内）中引入一个对y影响显著的；对集合中的变量进行检验，剔除影响变得不显著的；迭代式地进行引入和剔除，直到不能进行为止。选择衡量影响显著程度的统计量，通常用偏F统计量；适当选取引入变量的显著性水平in和剔除变量的out。引入新的变量后原来模型内影响显著的变量变得不显著，从而被剔除~自变量之间存在较强相关性的结果.MATLAB统计工具箱常用命令逐步回归

stepwise(x,y,inmodel,penter,premove)x~候选变量集合的n×k

数据矩阵（n是数据容量,k是变量数目）;y~因变量数据向量（n维）;Inmodel~初始模型中包括的候选变量集合的指标（矩阵x的列序数，缺省时设定为全部候选变量）;penter~引入变量的显著性水平（缺省时设定为0.05）;premove~剔除变量的显著性水平（缺省时设定为0.10）。输出交互式画面例教学评估为了考评教师的教学质量，教学研究部门设计了一个教学评估表。对学生对老师的课程进行打分。分值为1—5分（5分最好，1分最差）。x1—课程内容组织的合理性；x2—主要问题展开的合理性；x3—回答学生问题的有效性；x4—课下交流的有助性；x5—教科书的帮助性；x6—考试评分的公正性；y—对教师的总体评价。教师编号课程编号x1x2x3x4x5x6y12014.464.424.234.104.564.374.1122244.113.823.293.603.993.823.383301

3.583.313.243.764.393.753.17………………154244.244.384.354.484.154.504.33>>clearloaddata.txty=data(:,7);x=data(:,1:6);stepwise(x,y)例儿童的体重与身高和年龄序号体重(kg)身高(m)年龄序号体重(kg)身高(m)年龄127.11.348730.91.3910230.21.4910827.81.219324.01.146929.41.2610433.41.57111024.81.066524.91.1981136.51.6412624.31.1771229.11.449可能存在二次函数关系体重y身高x1体重y年龄x2例儿童的体重与身高和年龄初始结果最终结果例儿童的体重与身高和年龄初始结果最终结果[0，1]区间上的均匀随机数命令形式1：rand(N)功能：产生一个N*N的随机矩阵。引例抛一枚硬币10000次，如何模拟其正面的出现情况。命令形式2：rand(m,n)功能：产生一个m*n的随机矩阵。八、随机数正面x>=0.5反面x<0.5引例抛一枚硬币10000次，如何模拟其正面的出现情况。

fori=1:100

a(i)=sum(sum(round(rand(100))))/10000;enda

mx=max(a)

mn=min(a)ma=mean(a)白球x>0.7黑球x<=0.7例在箱子中有10个球，白球3个，黑球7个。P(白球)=0.3，P(黑球)=0.7。如何模拟该过程？

fori=1:100

a(i)=sum(sum(round(rand(100))))/10000;enda

mx=max(a)

mn=min(a)ma=mean(a)主成分分析思想九、主成分分析降维思想：高维到低维包含更多的、不重复的信息类型：总体主成分分析样本主成分分析1、总体主成分

定义设X1，X2，…，Xp

为某实际问题所涉及的p个随机变量。记X=(X1，X2，…,Xp)T，其协方差矩阵为

它是一个p阶非负定矩阵。设则有主成分析的目的第i个主成分：一般地，在约束条件

及下，求li使Var(Yi)达到最大，由此li所确定的称为X