Matlab学习系列17 数值计算概率篇.docx

Matlab学习系列17 数值计算概率篇.docx

《Matlab学习系列17 数值计算概率篇.docx》由会员分享，可在线阅读，更多相关《Matlab学习系列17 数值计算概率篇.docx（26页珍藏版）》请在冰豆网上搜索。

Matlab学习系列学习系列17数值计算概率篇数值计算概率篇17.数值计算概率篇一、计算组合数、排列数factorial（n）或prod（1:

n）nchoosek（n,k）factorial（n）/factorial（n-k）二、生成随机数1.rand（m,n）生成mn的服从0,1上均匀分布的随机数；用a+（b-a）.*rand（m,n）生成mn的服从a,b上均匀分布的随机数。

2.二项分布与正态分布随机数binornd（N,P,m,n）生成mn的服从二项分布B（N,P）的随机数；normrnd（MU,SIGMA,m,n）生成mn的服从正态分布N（MU,SIGMA2）的随机数；3.通用格式：

分布缩写+rnd（分布参数,m,n）或random（分布名或缩写,分布参数,m,n）可以用来生成mn该分布的随机数。

各种分布名见下图：

表1一维随机变量概率分布名称表分布缩写、分布名称函数说明beta或BetaBeta分布bino或Binomial二项分布chi2或Chisquare卡方分布exp或Exponential指数分布f或FF分布gam或GammaGAMMA分布geo或Geometric几何分布hyge或Hypergeometric超几何分布logn或Lognormal对数正态分布nbin或NegativeBinomial负二项式分布ncf或NoncentralF非中心F分布nct或Noncentralt非中心t分布ncx2或NoncentralChi-square非中心卡方分布norm或Normal正态分布poiss或Poisson泊松分布rayl或Rayleigh瑞利分布t或TT分布unif或Uniform均匀分布unid或DiscreteUniform离散均匀分布weib或Weibull威布尔分布4.使用randsample和randsrc函数生成指定离散分布随机数X=randsample（N,k,replace,w）N相当于1:

N,也可以是具有确定值的向量；k表示生成k个随机数；replace=true表示可重复，或false表示不可重复（默认）；w是权重向量。

X=randsrc（m,n,x;p）生成mn的随机矩阵，服从取值为向量x,对应概率为向量p的离散分布。

例1设离散型随机变量X服从如下分布：

X-2-1012P0.050.20.50.20.05生成服从35的该分布的随机数。

代码：

xvalue=-2-1012;xp=0.050.20.50.20.05;%调用randsample函数生成100个服从指定离散分布的随机数x=randsample（xvalue,15,true,xp）;reshape（x,35）%调用randsrc函数生成10*10的服从指定离散分布的随机数矩阵y=randsrc（3,5,xvalue;xp）运行结果：

ans=00100000-1-111001y=-1-111-1-10020-10-1005.已知概率密度函数，生成服从该分布的随机数例2设随机变量X的概率密度函数为（抛物线分布）：

调用crnd函数（来自MATLAB统计分析与应用40个案例分析作者：

谢中华），生成35个服从该分布的随机数。

代码：

pdffun=6*x*（1-x）;%密度函数表达式x=crnd（pdffun,01,3,5）运行结果：

x=0.31600.68660.27240.28160.12680.26810.84390.19480.79990.53830.73770.20400.49320.19480.69096.生成多元分布的随机数mrnd（N,P,m）多项分布，P为概率向量；mvnrnd（mu,sigma,m）多元正态分布，mu,sigma为n元向量；mvtrnd（C,df,m）多元t分布；wishrnd（sigma,df,m）Wishart分布；iwishrnd（sigma,df,m）逆Wishart分布；例3利用mvnrnd函数生成3组的二元正态分布随机数，其中分布的参数为代码：

mu=1020;sigma=13;316;xy=mvnrnd（mu,sigma,3）运行结果：

xy=11.833625.73859.034717.80269.603019.5821三、随机变量的概率密度函数及其图像概率密度函数，描述随机变量X在点x附近取值的可能性。

1.通用格式：

pdf（分布名或缩写,x,分布参数）返回该分布在X=x处的概率密度值；例如，Pk=pdf（bino,3,10,0.4）2.专用函数分布名缩写+pdf（x,分布参数）例如，binopdf（k,n,p）例4绘制卡方分布密度函数在自由度分别为1、5、15的图形。

代码：

x=0:

0.1:

30;y1=chi2pdf（x,1）;plot（x,y1,:

）holdony2=chi2pdf（x,3）;plot（x,y2,+）y3=chi2pdf（x,10）;plot（x,y3,o）axis（0,30,0,0.2）运行结果：

四、随机变量的分布函数分布函数定义为：

F（x）=PXx，表示随机变量X的取值落在（-,x）范围内的概率。

引入分布函数的目的，就是可以计算随机变量X的取值落在任意区间内的概率，例如，Pau=，故u=norminv（1-,0,1）其它分布的上分位数也是类似的。

N（0,1）的双侧分位数：

P|X|u/2=,故u/2=norminv（1-/2,0,1）注意：

-u/2=u1-/2（-u/2右侧的面积为1-/2），其它对称分布也成立，例如，t1-（n）=-t（n）.六、随机变量的数字特征注：

若要X中除NaN（非数）之外的数进行操作，加前缀nan,例如nanmean（X）.1.几种平均

（1）算术平均值（样本均值）适用于性质相同、单峰，且近似服从正态分布的定量数据；代码：

mean（X）（若X为矩阵，返回每列的均值）

（2）中位数代码：

median（X）（3）几何平均值适用于对比率数据的平均，并主要用于计算数据平均增长（变化）率；服从正偏态分布（较长右尾），特别是对数正态分布数据（取对数变换后服从正态分布）。

代码：

geomean（X）（4）调和平均值即先取倒数，再取算术平均值，再取倒数回来。

例如前半段时速60公里，后半段时速30公里（两段距离相等），则其平均速度为两者的调和平均数时速40公里。

在实际中，往往由于缺乏总体单位数的资料而不能直接计算算术平均数，这时需用调和平均法来求得平均数。

代码：

harmmean（X）（5）众数出现次数最多的数，代码：

mode（X）2.期望和方差

（1）样本数据mean（X）样本均值var（X）或var（X,0）样本方差var（X,1）方差std（X）或std（X,0）样本标准差Sstd（X,1）标准差

（2）常见分布的期望和方差通用格式：

E,D=分布名缩写+stat（分布参数）返回E为期望，D为方差；例如，E,D=normstat（MU,SIGMA）（3）已知离散型随机变量的分布律，求期望和方差例5设离散型随机变量X服从如下分布：

X-2-1012P0.30.10.20.10.3求E（X）,E（X2-1）,D（X）.代码：

X=-2-1012;p=0.30.10.20.10.3;EX=sum（X.*p）%或EX=X*pY=X.2-1;EY=sum（Y.*p）XX=X.2;EXX=sum（XX.*p）;DX=EXX-EX2运行结果：

EX=0EY=1.6000DX=2.60003.极差、偏度、峰度极差，数据的最大值与最小值之差。

代码：

range（X）偏度，是数据关于均值不对称的指标。

若偏度为负，说明均值左边的数据比右边的数据更散；若偏度为正，说明均值右边的数据比左边的数据更散，正态分布的偏度为0.代码：

skewness（X）峰度，是用来反映数据的分布曲线顶端尖峭或扁平程度的指标。

代码：

kurtosis（X）4.矩k阶中心矩m=moment（X,k）k阶原点矩可以使用下面自编的OriginMoment.mfunctiony=OriginMoment（X,k）;%X为样本矩阵n,m=size（X）;y=zeros（1,m）;forii=1:

m;y（ii）=sum（X（:

ii）.k）/n;end5.协方差（矩阵）、相关系数（矩阵）对于样本数据矩阵：

Matlab将矩阵X的每一列Xj,j=1,n作为一个随机变量的样本，每一行xj1,xj2,xjn,j=1,2,.,m作为n个随机变量的联合分布的一个样本。

由于矩阵X给出的只是随机变量的样本数据，并不知道这些随机变量的（联合）概率分布，因此是不能计算出这些随机变量的总体期望、方差或协方差的，而只能计算出它们的一个无偏估计，即样本均值、样本方差与样本协方差。

样本协方差公式：

其中称为随机变量X1,Xn的协方差矩阵（实对称、非负定）。

协方差矩阵是一维随机变量方差、二维随机向量协方差向高维随机向量的推广，常应用于在主成分分析中。

相关系数（矩阵）是协方差（矩阵）的标准化，反映了随机变量两两之间的线性关系的强弱程度。

代码：

cov（X）返回样本矩阵X的协方差矩阵；corrcoef（X）返回样本矩阵X的相关系数矩阵；例6随机生成样本数据矩阵X，计算X的协方差矩阵、相关系数矩阵。

代码：

M=5;N=3;X=10.*rand（M,N）CovX=cov（X）Cov12=cov（X（:

1）,X（:

2）Cov11=cov（X（:

1）,X（:

1）var（X（:

1）RhoX=corrcoef（X）运行结果：

X=1.06658.68694.31419.61900.84449.10650.04633.99781.81857.74912.59872.63808.17308.00071.4554CovX=19.6061-6.38125.3900-6.381211.6214-5.58945.3900-5.58949.7937Cov12=19.6061-6.3812-6.381211.6214Cov11=19.606119.606119.606119.6061VarX1=19.6061RhoX=1.0000-0.42270.3890-0.42271.0000-0.52390.3890-0.52391.0000七、频率直方图1.统计数值或字符出现的频数、频率设X为数值型或字符型向量，代码：

tabulate（X）;输出X的值、频数、频率的三列表。

注意：

若X中的值都是非负整数，则min（X）与max（X）之间的任一整数都统计，未出现的频数=0.2.经验分布函数实际上是累积频率直方图曲线。

依据样本以频率估计概率，当n充分大时，是总体分布函数的近似。

例7生成指数分布随机数，计算其经验分布函数值并绘图，并与其分布函数做对比。

代码：

X=exprnd（3,100,1）;fp,xp=ecdf（X）;%计算经验分布函数值subplot（1,2,1）ecdfhist（fp,xp,20）;%绘制频率直方图subplot（1,2,2）plot（xp,fp,r）%绘制经验分布函数图形（连线型）holdonh,stats=cdfplot（X）%绘制经验分布函数图形（阶梯型）x=0:

0.01:

15;y=expcdf（x,3）;plot（x,y,g）;legend（连线型经验分布函数,阶梯型经验分布函数,分布函数）;运行结果：

h=177.0125（图形句柄）stats=min:

0.0180max:

27.9801mean:

3.2968median:

2.1667std:

3.29690510152000.050.10.150.20.250.30.350510152000.10.20.30.40.50.60.70.80.91xF（x）EmpiricalCDF连线型经验分布函数阶梯型经验分布函数分布函数3.直方图定量数据常用直方图来展示某变量取值的分布，利用直方图可以估计总体的概率密度。

例8随机生成服从F（3,5）分布的样本数据，绘制直方图并与该分布的真实概率密度曲线做对比。

代码：

n=5000;X=frnd（3,5,n,1）;%生成F分布随机数m=150;%分组区间数a=min（X）;b=max（X）;d=（b-a）/m;%分组宽度r,xout=hist（X,a:

d:

b）;%或者r,xout=hist（X,m）;%计算直方图数据,r返回每段的频数,xout返回每段的中心位置f=r./（n*d）;%计算频率bar（xout,f）;%绘制频率直方图holdonx=0:

0.01:

10;y=fpdf（x,3,5）;plot（x,y,k-）;axis（01001）;title（频率密度直方图）;运行结果：

注：

也可以用例7中方法绘制频率直方图：

fp,xp=ecdf（X）;%计算经验分布函数值ecdfhist（fp,xp,20）;%绘制频率直方图4.绘制分布的概率图形normplot（X）绘制X的正态分布概率图形（若X为矩阵，则针对每列绘制）；可以用来验证数据X是否服从正态分布，若是，则近似一条直线。

类似的还有，weibplot（X）绘制X的威布尔分布概率图形。

例9随机生成服从正态分布的随机数据，绘制其正态概率分布图形，验证数据的正态性。

代码：

X=normrnd（5,1.44,100,1）;normplot（X）运行结果：

024681000.10.20.30.40.50.60.70.80.91频率密度直方图八、箱线图又称为盒形图。

在一条数轴上，以数据的上下四分位数（Q1-Q3）为界画一个矩形盒子（中间50%的数据落在盒内）；在数据的中位数位置画一条线段为中位线；默认延长线不超过盒长的1.5倍，之外的点认为是异常值（用+标记）。

盒形图的主要应用就是，剔除数据的异常值、判断数据的偏态和尾重。

代码：

boxplot（X,notch,sym,vert,whis）当notch=1时，产生一凹盒图，notch=0时产生一矩箱图；sym表示异常值的标记符号，默认值为“+”；当vert=0时，生成水平盒图，vert=1时，生成竖直盒图（默认值vert=1）；whis定义“须线”的长度，默认值为1.5.例10绘制箱线图。

代码：

x1=normrnd（5,1,100,1）;x2=normrnd（6,1,100,1）;1234567890.0030.010.020.050.100.250.500.750.900.950.980.990.997DataProbabilityNormalProbabilityPlotx=x1x2;boxplot（x,0,g+,1,1.5）运行结果：

注：

若有异常值，应先去掉异常值，再重新绘制修正的箱线图。

九、参数估计刻画总体某方面概率特性的量，称为参数；当某参数未知时，从总体抽取样本，用某种方法（矩估计、最大似然估计）对该未知参数进行估计，就是参数估计。

通过构造样本的函数，给出未知参数的估计值（点估计）、取值范围（区间估计）。

1.通过格式（基于最大似然估计）phat,pci=分布名缩写+fit（X,alpha）phat返回参数的点估计值，pci返回区间估计的两端点；alpha为显著水平，即估计结果具有（1-alpha）100%的置信度。

注意：

对于二项分布的参数估计需要带参数“N”表示试验次数：

phat,pci=binofit（x,N,alpha）2.或者使用phat,pci=mle（分布名缩写,X,alpha）345678912对于二项分布为：

phat,pci=mle（分布名缩写,X,alpha,N）3.求最大对数似然函数值logL,info=分布名缩写+like（给定参数,X）其中，给定参数为行向量形式，输入的是参数的极大似然估计值；info返回回Fisher逆信息矩阵，其对角线为为相应参数的渐近方差。

例11某厂生产的滚珠随机抽取10个，测得滚珠直径（mm）如下：

15.1414.8115.1115.2615.0815.1715.1214.9515.0514.87设测定值总体服从N（,2）,和未知，分别求和的最大似然估计以及95%的置信区间。

代码：

X=15.1414.8115.1115.2615.0815.1715.1214.9515.0514.87;muhat,sigmaheat,muci,sigmaci=normfit（X,0.05）运行结果：

muhat=15.0560sigmahat=0.1397muci=14.956115.1559sigmaci=0.09610.2550logL=-5.9933info=0.00200.00000.00000.0011十、假设检验实际中，我们只能抽取部分样本，做统计分析得到统计结果，进一步推断总体的特征，但是这种推断必然有可能犯错，犯错的概率为多少时应该接受这种推断呢？

为此，统计学家就开发了一些统计方法进行统计检定，通过把所得到的统计检定值，与统计学家树立了一些随机变量的概率分布进行对比，我们可以知道在百分之多少的机遇下会得到目前的结果。

倘若经比较后发现，涌现这结果的机率很少，即是说，是在时机很少、很罕有的情况下才出现；那我们便可以有信念地说，这不是巧合，该推断结果是具有统计学上的意义的。

否则，就是推断结果不具有统计学意义。

假设检验是基于反证法思想：

先提出原假设（H0），再用适当的统计方法确定假设成立的可能性（P值）大小，如可能性小（P），则认为原假设不成立，若可能性大，则还不能认为备择假设（H1）成立。

原假设与备择假设是是完备且相互独立的事件组，一般，原假设（H0）研究者想收集证据予以反对的假设；备择假设（H1）研究者想收集证据予以支持的假设；假设检验的P值是由检验统计量的样本观察值得出的原假设可被拒绝的最小显著水平。

（1）双侧检验（H1:

0）I.原假设H0:

=0,备择假设H1:

0;.根据样本数据计算出统计量t的观察值t0;.P值=P|t|t0|=t0的双侧尾部的面积；.若P值（在右尾部分），则在显著水平下拒绝H0;若P值，则在显著水平下接受H0;注意：

为临界值，看P值在不在阴影部分（拒绝域），空白部分为接受域。

（2）左侧检验（H1:

0）I.原假设H0:

0,备择假设H1:

0;.根据样本数据计算出统计量t的观察值t0（，则在显著水平下接受H0;（3）右侧检验（H1:

0）I.原假设H0:

0,备择假设H1:

0;.根据样本数据计算出统计量t的观察值t0（0）;.P值=Ptt0=t0的右侧尾部的面积；.若P值（在右尾部分），则在显著水平下拒绝H0;若P值，则在显著水平下接受H0;1.总体标准差已知，单个正态总体均值的U检验h,p,muci,zval=ztest（X,mu,sigma,alpha,tail）p返回p值，h=1或palpha时，拒绝H0;h=0或palpha时接受H0.muci返回总体均值的1-alpha置信区间；zval返回检验统计量的观测值；tail=both或right或left,用来规定双侧、右侧、左侧检验，默认是both.例12某切割机正常工作时，切割的金属棒的长度服从N（100,4）,随机抽取15根，测得长度（mm）如下：

9710210511299103102941009510598102100103假设总体方差不变，检验该切割机是否正常，即总体均值是否等于100mm？

显著性水平取=0.05.H0:

=0=100,H1:

0代码：

X=9710210511299103102941009510598102100103;h,p,muci,zval=ztest（X,100,2,0.05）%h=1或p0.05时，拒绝原假设；h=0或p0.05时，接受原假设h1,p1,muci1,zval1=ztest（X,100,2,0.05,right）运行结果：

h=1p=0.0282muci=100.1212102.1455zval=2.1947（U统计量值）h1=1p1=0.0141uci1=100.2839Infzval1=2.1947由于p值=0.02820经检验，p值=0.0141=0.05,拒绝原假设H0,故结论是：

在显著水平=0.05下，总体均值大于100.2.t检验t检验，用来检验“单样本均值与已知均值、两独立样本均值、配对设计资料的均值”是否存在差异，这种差异是否能推论至总体。

t检验适用于样本含量较小（比如n=0.05,接受原假设，故结论是：

在显著水平=0.05下，每包化肥的平均质量等于50kg.例14甲、乙两台机床加工同一种产品，从这两台机床加工的产品中随机抽取若干件，测得产品直径（mm）如下：

甲机床：

20.120.019.320.620.219.920.019.919.119.9乙机床：

18.619.120.020.020.019.719.919.620.2设甲、乙两机床加工的产品直径分别服从N（1,12）和N（2,22），试比较两台机床加工的产品直径是否有差异？

显著性水平取=0.05.H0:

1=2,H1:

12代码：

X=20.120.019.320.620.219.920.019.919.119.9;Y=18.619.120.020.020.019.719.919.620.2;h,p,muci,stats=ttest2（X,Y,0.05,both,equal）运行结果：

h=0p=0.3191muci=-0.23460.6791stats=tstat:

1.0263df:

17sd:

0.4713经检验，p值=0.3191=0.05,接受原假设，故结论是：

在显著水平=0.05下，两台机床加工的产品直径没有显著差异。

3.总体均值未知，单个正态总体方差的2检验卡方的统计原理，是取观察频数与期望频数相比较。

当观察频数与期望频数完全一致时，值为0；观察频数与期望频数越接近，两者之间的差异越小，值越小；观察频数与期望频数差别越大，两者之间的差异越大，值越大。

一旦值大于某一个临界值，即可获得显著的统计结论。

配对资料2检验研究两变量之间的关联性或差别性问题。

如果两变量无关联（即相互独立），说明对于其中一个变量而言，另一变量多项分类次数上的变化是在无差范围之内；如果两变量有关联，说明二者之间有交互作用存在。

两变量之间是否相互关联，也可以看作差别分析，例如说性别与成绩有关，其实意思就是说不同性别的人他们的成绩有差别。

这里用2检验来检验方差与某方差值是否有差别，代码：

h,p,muci,stats=vartest（X,var0,alpha,tail）例15还是例13中的样本数据，检验每包化肥质量的方差是否等于1.5？

显著性水平取=0.05.H0:

2=02=1.5,H1:

202代码：

X=49.450.550.751.749.847.949.2