1、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. 通用格式：分布缩写+
2、rnd(分布参数, m,n)或random(分布名或缩写, 分布参数, m,n)可以用来生成mn该分布的随机数。各种分布名见下图：表1 一维随机变量概率分布名称表分布缩写、分布名称函数说明beta 或 Beta Beta分布bino 或 Binomial 二项分布chi2 或 Chisquare 卡方分布exp 或 Exponential 指数分布f 或 F F分布gam 或 Gamma GAMMA分布geo 或 Geometric 几何分布hyge 或 Hypergeometric 超几何分布logn 或 Lognormal 对数正态分布nbin 或 Negative Binomial 负二
3、项式分布ncf 或 Noncentral F 非中心F分布nct 或 Noncentral t 非中心t分布ncx2 或 Noncentral Chi-square非中心卡方分布norm 或 Normal 正态分布poiss 或 Poisson 泊松分布rayl 或 Rayleigh 瑞利分布t 或 T T分布unif 或 Uniform 均匀分布unid 或 Discrete Uniform 离散均匀分布weib 或 Weibull 威布尔分布4. 使用randsample和randsrc函数生成指定离散分布随机数X=randsample(N, k, replace, w)N相当于1:N,
4、也可以是具有确定值的向量；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 -1 0 1 2; xp = 0.05 0.2 0.5 0.2 0.05; % 调用randsample函数生成100个服从指定离散分布的随机数x = randsample(xvalue, 15, true,
5、xp);reshape(x,3 5)% 调用randsrc函数生成10*10的服从指定离散分布的随机数矩阵y = randsrc(3,5,xvalue;xp)运行结果：ans = 0 0 1 0 00 0 0 -1 -11 1 0 0 1y = -1 -1 1 1 -1-1 0 0 2 0-1 0 -1 0 05. 已知概率密度函数，生成服从该分布的随机数例2 设随机变量X的概率密度函数为（抛物线分布）：调用crnd函数（来自MATLAB统计分析与应用 40个案例分析作者：谢中华），生成35个服从该分布的随机数。代码：pdffun = 6*x*(1-x); % 密度函数表达式x = crnd(
6、pdffun,0 1,3,5)运行结果：x = 0.3160 0.6866 0.2724 0.2816 0.12680.2681 0.8439 0.1948 0.7999 0.53830.7377 0.2040 0.4932 0.1948 0.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函
7、数生成3组的二元正态分布随机数，其中分布的参数为代码：mu = 10 20; sigma = 1 3; 3 16; xy = mvnrnd(mu, sigma, 3)运行结果：xy = 11.8336 25.73859.0347 17.80269.6030 19.5821三、随机变量的概率密度函数及其图像概率密度函数，描述随机变量X在点x附近取值的可能性。1. 通用格式：pdf(分布名或缩写, x, 分布参数)返回该分布在X=x处的概率密度值；例如，Pk=pdf(bino,3, 10, 0.4)2. 专用函数分布名缩写+pdf(x, 分布参数)例如，binopdf (k, n, p)例4绘制卡
8、方分布密度函数在自由度分别为1、5、15的图形。代码：x=0:0.1:30;y1=chi2pdf(x,1); plot(x,y1,:)hold ony2=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/
9、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) 几何平均值适用于对比率数据的平均，并主要用于计算数据平均增长（变化）率；服从正偏态分布（较长右尾），特别是对数正态分布数据（取对数变
10、换后服从正态分布）。代码：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,
11、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 -1 0 1 2;p=0.3 0.1 0.2 0.1 0.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. 极差、偏度、峰度
12、极差，数据的最大值与最小值之差。代码：range(X)偏度，是数据关于均值不对称的指标。若偏度为负，说明均值左边的数据比右边的数据更散；若偏度为正，说明均值右边的数据比左边的数据更散，正态分布的偏度为0. 代码：skewness(X)峰度，是用来反映数据的分布曲线顶端尖峭或扁平程度的指标。代码：kurtosis(X)4. 矩k阶中心矩m=moment(X, k)k阶原点矩可以使用下面自编的OriginMoment.mfunction y=OriginMoment(X,k);% X为样本矩阵 n,m=size(X); y=zeros(1,m);for ii=1:m;y(ii)=sum(X(:,i
13、i).k)/n;end5. 协方差（矩阵）、相关系数（矩阵）对于样本数据矩阵：Matlab将矩阵X的每一列Xj, j=1,n作为一个随机变量的样本，每一行xj1,xj2, , xjn, j=1,2,.,m作为n个随机变量的联合分布的一个样本。由于矩阵X给出的只是随机变量的样本数据，并不知道这些随机变量的（联合）概率分布，因此是不能计算出这些随机变量的总体期望、方差或协方差的，而只能计算出它们的一个无偏估计，即样本均值、样本方差与样本协方差。样本协方差公式：, 其中称为随机变量X1, , Xn的协方差矩阵（实对称、非负定）。协方差矩阵是一维随机变量方差、二维随机向量协方差向高维随机向量的推广，常
14、应用于在主成分分析中。相关系数（矩阵）是协方差（矩阵）的标准化，反映了随机变量两两之间的线性关系的强弱程度。代码：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.0665 8.6869 4.3141 9.6190 0.844
15、4 9.1065 0.0463 3.9978 1.8185 7.7491 2.5987 2.6380 8.1730 8.0007 1.4554CovX = 19.6061 -6.3812 5.3900-6.3812 11.6214 -5.58945.3900 -5.5894 9.7937Cov12 = 19.6061 -6.3812-6.3812 11.6214Cov11 = 19.6061 19.606119.6061 19.6061VarX1 = 19.6061RhoX = 1.0000 -0.4227 0.3890-0.4227 1.0000 -0.52390.3890 -0.5239
16、1.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); %
17、 绘制频率直方图subplot(1,2,2)plot(xp,fp,r) % 绘制经验分布函数图形(连线型）hold onh,stats=cdfplot(X) % 绘制经验分布函数图形(阶梯型）x=0:0.01:15;y=expcdf(x,3);plot(x,y,g);legend(连线型经验分布函数,阶梯型经验分布函数,分布函数);运行结果：h = 177.0125（图形句柄）stats = min: 0.0180 max: 27.9801mean: 3.2968 median: 2.1667 std: 3.29693. 直方图定量数据常用直方图来展示某变量取值的分布，利用直方图可以估计总体的
18、概率密度。例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); % 绘制频率直方图hold onx=0:0.01:10;y=fpdf(x,3,5);plot(x,y,k-)
19、;axis(0 10 0 1);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
20、(X)运行结果：八、箱线图又称为盒形图。在一条数轴上，以数据的上下四分位数（Q1-Q3）为界画一个矩形盒子（中间50%的数据落在盒内）；在数据的中位数位置画一条线段为中位线；默认延长线不超过盒长的1.5倍，之外的点认为是异常值（用+标记）。盒形图的主要应用就是，剔除数据的异常值、判断数据的偏态和尾重。代码：boxplot(X,notch,sym,vert,whis)当notch=1时，产生一凹盒图，notch=0时产生一矩箱图；sym表示异常值的标记符号，默认值为“+”；当vert=0时，生成水平盒图，vert=1时，生成竖直盒图（默认值vert=1）；whis定义“须线”的长度，默认值为1.
21、5.例10 绘制箱线图。代码：x1 = normrnd(5,1,100,1);x2 = normrnd(6,1,100,1);x = x1 x2;boxplot(x,0,g+,1,1.5)运行结果：注：若有异常值，应先去掉异常值，再重新绘制修正的箱线图。九、参数估计刻画总体某方面概率特性的量，称为参数；当某参数未知时，从总体抽取样本，用某种方法（矩估计、最大似然估计）对该未知参数进行估计，就是参数估计。通过构造样本的函数，给出未知参数的估计值（点估计）、取值范围（区间估计）。1. 通过格式（基于最大似然估计）phat,pci=分布名缩写+fit(X, alpha)phat返回参数的点估计值，p
22、ci返回区间估计的两端点；alpha为显著水平，即估计结果具有(1-alpha)100%的置信度。注意：对于二项分布的参数估计需要带参数“N”表示试验次数：phat, pci= binofit (x, N, alpha)2. 或者使用phat,pci=mle(分布名缩写, X, alpha)对于二项分布为：phat,pci=mle(分布名缩写, X, alpha, N)3. 求最大对数似然函数值logL, info=分布名缩写+like(给定参数, X)其中，给定参数为行向量形式，输入的是参数的极大似然估计值；info返回回Fisher逆信息矩阵，其对角线为为相应参数的渐近方差。例11某厂生产
23、的滚珠随机抽取10个，测得滚珠直径(mm)如下：15.14 14.81 15.11 15.26 15.0815.17 15.12 14.95 15.05 14.87设测定值总体服从N(, 2), 和未知，分别求和的最大似然估计以及95%的置信区间。代码：X=15.14 14.81 15.11 15.26 15.08 15.17 15.12 14.95 15.05 14.87;muhat,sigmaheat,muci,sigmaci=normfit(X, 0.05)运行结果：muhat = 15.0560sigmahat = 0.1397muci = 14.956115.1559sigmaci
24、= 0.09610.2550logL = -5.9933info = 0.0020 0.0000 0.0000 0.0011十、假设检验实际中，我们只能抽取部分样本，做统计分析得到统计结果，进一步推断总体的特征，但是这种推断必然有可能犯错，犯错的概率为多少时应该接受这种推断呢？为此，统计学家就开发了一些统计方法进行统计检定，通过把所得到的统计检定值，与统计学家树立了一些随机变量的概率分布进行对比，我们可以知道在百分之多少的机遇下会得到目前的结果。倘若经比较后发现，涌现这结果的机率很少，即是说，是在时机很少、很罕有的情况下才出现；那我们便可以有信念地说，这不是巧合，该推断结果是具有统计学上的意义
25、的。否则，就是推断结果不具有统计学意义。假设检验是基于反证法思想：先提出原假设（H0），再用适当的统计方法确定假设成立的可能性（P值）大小，如可能性小（P），则认为原假设不成立，若可能性大，则还不能认为备择假设（H1）成立。原假设与备择假设是是完备且相互独立的事件组，一般，原假设（H0）研究者想收集证据予以反对的假设；备择假设（H1）研究者想收集证据予以支持的假设；假设检验的P值是由检验统计量的样本观察值得出的原假设可被拒绝的最小显著水平。（1）双侧检验（H1: 0）I. 原假设H0: =0, 备择假设H1: 0;. 根据样本数据计算出统计量t的观察值t0;. P值 = P|t| |t0| =
26、 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值 = Pt t0 = t0的右侧尾部的面积；. 若P值（在右尾部分），则在显著水平下拒绝H0;若P值，则在显著水平下接受H0;1. 总体标准差已知，单
27、个正态总体均值的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)如下：97 102 105 112 99 103 102 94 100 95 105 98 102 100 103假设总体方差不变，检验该切割机是否正常，即总体均值是否等于100mm？显著性水平取=0.05.H0: =0=100, H1: 0代码：X=97 102 105 112 99 103 102 94 100 95 105 98 102 100 103;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.1212 10
