1、17. 数值计算概率篇一、计算组合数、排列数factorial(n) 或 prod(1:n)!nnchoosek(n,k)kCfactorial(n)/factorial(n-k)nA二、生成随机数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. 通用格式：分布缩写+rn
2、d(分布参数, 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 或 Negat
3、ive Binomial 负二项式分布 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=r
4、andsample(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 -1 0 1 2P 0.05 0.2 0.5 0.2 0.05生成服从 35 的该分布的随机数。代码：xvalue = -2 -1 0 1 2; xp = 0.05 0.2 0.5 0.2 0.05; % 调用
5、randsample 函数生成 100 个服从指定离散分布的随机数x = randsample(xvalue, 15, true, 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 的概率密度函数为（抛物线分布）： 6(1), 01() xxf其 他调用 c
6、rnd 函数（来自MATLAB 统计分析与应用 40 个案例分析作者：谢中华） ，生成 35 个服从该分布的随机数。代码：pdffun = 6*x*(1-x); % 密度函数表达式x = crnd(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
7、 为 n 元向量；mvtrnd(C, df, m)多元 t 分布；wishrnd(sigma,df,m)Wishart 分布；iwishrnd(sigma,df, m)逆 Wishart 分布；例 3 利用 mvnrnd 函数生成 3 组的二元正态分布随机数，其中分布的参数为 10=2316，代码：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 附近取
8、值的可能性。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,:)hold ony2=chi2pdf(x,3);plot(x,y2,+)y3=chi2pdf(x,10);plot(x,y3,o)axis(0,30,0,0.2)运行结果：0 5 10 15 20
9、25 3000.020.040.060.080.10.120.140.160.180.2四、随机变量的分布函数分布函数定义为：F(x)=PX x，表示随机变量 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） ，其它对称分布也成立，例如，t 1-(n)=- t(n).六、随机变量的数字特征
10、注：若要 X 中除 NaN（非数）之外的数进行操作，加前缀 nan, 例如 nanmean(X).1. 几种平均(1) 算术平均值（样本均值）适用于性质相同、单峰，且近似服从正态分布的定量数据；代码：mean(X) （若 X 为矩阵，返回每列的均值）(2) 中位数代码：median(X)(3) 几何平均值 1niGx适用于对比率数据的平均，并主要用于计算数据平均增长（变化）率；服从正偏态分布（较长右尾） ，特别是对数正态分布数据（取对数变换后服从正态分布） 。代码：geomean(X)(4) 调和平均值 1niHx即先取倒数，再取算术平均值，再取倒数回来。例如前半段时速 60公里，后半段时速
11、30 公里（两段距离相等） ，则其平均速度为两者的调和平均数时速 40 公里。在实际中，往往由于缺乏总体单位数的资料而不能直接计算算术平均数，这时需用调和平均法来求得平均数。代码：harmmean(X)(5) 众数出现次数最多的数，代码：mode(X)2. 期望和方差(1) 样本数据mean(X)样本均值 1nixvar(X)或 var(X,0)样本方差 221=()niiSxvar(X,1)方差 21()niiDxstd(X)或 std(X,0)样本标准差 Sstd(X,1)标准差 (2) 常见分布的期望和方差通用格式：E,D=分布名缩写+stat(分布参数)返回 E 为期望，D 为方差；例
12、如，E,D=normstat(MU,SIGMA)(3) 已知离散型随机变量的分布律，求期望和方差例 5 设离散型随机变量 X 服从如下分布：X -2 -1 0 1 2P 0.3 0.1 0.2 0.1 0.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 = 0 EY = 1.6000 DX = 2.60003. 极差、偏度、峰度极差，数据的最大值与最小值
13、之差。代码：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(:,ii).k)/n;en
14、d5. 协方差（矩阵） 、相关系数（矩阵）对于样本数据矩阵：12121212n nmmnxxXX Matlab 将矩阵 X 的每一列 Xj, j=1,n 作为一个随机变量的样本，每一行 xj1, xj2, , xjn, j=1,2,.,m 作为 n 个随机变量的联合分布的一个样本。由于矩阵 X 给出的只是随机变量的样本数据，并不知道这些随机变量的（联合）概率分布，因此是不能计算出这些随机变量的总体期望、方差或协方差的，而只能计算出它们的一个无偏估计，即样本均值、样本方差与样本协方差。样本协方差公式： 1(,)()niiiCOVXYxy, 其中12212nnncc (,)ijijcCOVX称为随
15、机变量 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),
16、X(:,1)var(X(:,1)RhoX = corrcoef(X)运行结果：X = 1.0665 8.6869 4.31419.6190 0.8444 9.10650.0463 3.9978 1.81857.7491 2.5987 2.63808.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 =
17、19.6061RhoX = 1.0000 -0.4227 0.3890-0.4227 1.0000 -0.52390.3890 -0.5239 1.0000七、频率直方图1. 统计数值或字符出现的频数、频率设 X 为数值型或字符型向量，代码：tabulate(X); 输出 X 的值、频数、频率的三列表。注意：若 X 中的值都是非负整数，则 min(X)与 max(X)之间的任一整数都统计，未出现的频数=0.2. 经验分布函数 1, ()nnxxF 中 的 个 数实际上是累积频率直方图曲线。依据样本以频率估计概率，当 n 充分大时， 是总体分布函数 的近似。()nx()x例 7 生成指数分布随机
18、数，计算其经验分布函数值并绘图，并与其分布函数做对比。代码：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) % 绘制经验分布函数图形( 连线型）hold onh,stats=cdfplot(X) % 绘制经验分布函数图形(阶梯型）x=0:0.01:15;y=expcdf(x,3);plot(x,y,g);legend(连线型经验分布函数 ,阶梯型经验分布函数, 分布函数);运行结果：h = 177.0125 （
19、图形句柄）stats = min: 0.0180 max: 27.9801mean: 3.2968 median: 2.1667 std: 3.29690 5 10 15 2000.050.10.150.20.250.30.350 5 10 15 2000.10.20.30.40.50.60.70.80.91xF(x)Empirical CDF函函函函函函函函函函函函函函函函函函函函函函3. 直方图定量数据常用直方图来展示某变量取值的分布，利用直方图可以估计总体的概率密度。例 8 随机生成服从 F(3,5)分布的样本数据，绘制直方图并与该分布的真实概率密度曲线做对比。代码：n=5000;X=f
20、rnd(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-);axis(0 10 0 1);title(频率密度直方图);运行结果：0 2 4 6 8 1000.10.20.30
21、.40.50.60.70.80.91 函函函函函函函注：也可以用例 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);normp
22、lot(X)运行结果：1 2 3 4 5 6 7 8 90.0030.010.020.050.100.250.500.750.900.950.980.990.997DataProbabilityNormal Probability Plot八、箱线图又称为盒形图。在一条数轴上，以数据的上下四分位数（Q 1-Q3）为界画一个矩形盒子（中间 50%的数据落在盒内） ；在数据的中位数位置画一条线段为中位线；默认延长线不超过盒长的 1.5 倍，之外的点认为是异常值（用+标记） 。盒形图的主要应用就是，剔除数据的异常值、判断数据的偏态和尾重。代码：boxplot(X,notch,sym,vert,whi
23、s)当 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);x = x1 x2;boxplot(x,0,g+,1,1.5)运行结果：34567891 2注：若有异常值，应先去掉异常值，再重新绘制修正的箱线图。九、参数估计刻画总体某方面概率特性的量，称为参数；当某参数未知
24、时，从总体抽取样本，用某种方法（矩估计、最大似然估计）对该未知参数进行估计，就是参数估计。通过构造样本的函数，给出未知参数的估计值（点估计） 、取值范围（区间估计） 。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)对于二项分
25、布为：phat,pci=mle(分布名缩写 , X, alpha, N)3. 求最大对数似然函数值logL, info=分布名缩写+like( 给定参数 , X)其中，给定参数为行向量形式，输入的是参数的极大似然估计值；info 返回回 Fisher 逆信息矩阵，其对角线为为相应参数的渐近方差。例 11 某厂生产的滚珠随机抽取 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 1
26、4.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.0560 sigmahat = 0.1397muci = 14.956115.1559sigmaci = 0.09610.2550logL = -5.9933info = 0.0020 0.00000.0000 0.0011十、假设检验实际中，我们只能抽取部分样本，做统计分析得到统计结果，进一步推断总体的特征，但是这种推断必然有可能犯错，犯错的概率为多少时应该接受这种推
27、断呢？为此，统计学家就开发了一些统计方法进行统计检定，通过把所得到的统计检定值，与统计学家树立了一些随机变量的概率分布进行对比，我们可以知道在百分之多少的机遇下会得到目前的结果。倘若经比较后发现，涌现这结果的机率很少，即是说，是在时机很少、很罕有的情况下才出现；那我们便可以有信念地说，这不是巧合，该推断结果是具有统计学上的意义的。否则，就是推断结果不具有统计学意义。假设检验是基于反证法思想：先提出原假设（H 0） ，再用适当的统计方法确定假设成立的可能性（P 值）大小，如可能性小（P ） ，则认为原假设不成立，若可能性大，则还不能认为备择假设（H 1）成立。原假设与备择假设是是完备且相互独立的
28、事件组，一般，原假设（H 0）研究者想收集证据予以反对的假设；备择假设（H 1）研究者想收集证据予以支持的假设；假设检验的 P 值是由检验统计量的样本观察值得出的原假设可被拒绝的最小显著水平。（1）双侧检验（H 1: 0）I. 原假设 H0: =0, 备择假设 H1: 0;. 根据样本数据计算出统计量 t 的观察值 t0;. P 值 = P|t| |t0| = t0 的双侧尾部的面积；. 若 P 值 （在右尾部分） ，则在显著水平 下拒绝 H0;若 P 值 ，则在显著水平 下接受 H0;注意： 为临界值，看 P 值在不在阴影部分（拒绝域） ，空白部分为接受域。（2）左侧检验（H 1: ，则在显
29、著水平 下接受 H0;（3）右侧检验（H 1: 0）I. 原假设 H0: 0, 备择假设 H1: 0;. 根据样本数据计算出统计量 t 的观察值 t0（ 0）;. P 值 = Pt t0 = 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
30、置信区间；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
31、 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 102.1455zval = 2.1947 （U 统计量值）h1 = 1p1 = 0.0141uci1 = 100.2839 Infzval1 = 2.1947由于 p 值=0.02820经检验，p 值=0.0141=0.05, 接受原假设，故结论是：在显著水平=0.05 下，每包化
32、肥的平均质量等于 50kg. 例 14 甲、乙两台机床加工同一种产品，从这两台机床加工的产品中随机抽取若干件，测得产品直径(mm)如下：甲机床：20.1 20.0 19.3 20.6 20.2 19.9 20.0 19.9 19.1 19.9乙机床：18.6 19.1 20.0 20.0 20.0 19.7 19.9 19.6 20.2设甲、乙两机床加工的产品直径分别服从 N(1, 12)和 N(2, 22)，试比较两台机床加工的产品直径是否有差异？显著性水平取 =0.05.H0: 1=2, H1: 1 2代码：X=20.1 20.0 19.3 20.6 20.2 19.9 20.0 19.9
33、 19.1 19.9;Y=18.6 19.1 20.0 20.0 20.0 19.7 19.9 19.6 20.2;h,p,muci,stats=ttest2(X,Y,0.05,both,equal)运行结果：h = 0p = 0.3191muci = -0.2346 0.6791stats = tstat: 1.0263df: 17sd: 0.4713经检验，p 值=0.3191=0.05, 接受原假设，故结论是：在显著水平=0.05 下，两台机床加工的产品直径没有显著差异。3. 总体均值未知，单个正态总体方差的 2 检验卡方的统计原理，是取观察频数与期望频数相比较。当观察频数与期望频数完全
34、一致时， 值为 0；观察频数与期望频数越接近，2两者之间的差异越小， 值越小；观察频数与期望频数差别越大，两者之间的差异越大， 值越大。一旦 值大于某一个临界值，即22可获得显著的统计结论。配对资料 2 检验研究两变量之间的关联性或差别性问题。如果两变量无关联（即相互独立） ，说明对于其中一个变量而言，另一变量多项分类次数上的变化是在无差范围之内；如果两变量有关联，说明二者之间有交互作用存在。两变量之间是否相互关联，也可以看作差别分析，例如说性别与成绩有关，其实意思就是说不同性别的人他们的成绩有差别。这里用 2 检验来检验方差与某方差值是否有差别，代码：h,p,muci,stats=varte
35、st(X,var0, alpha, tail)例 15 还是例 13 中的样本数据，检验每包化肥质量的方差是否等于1.5？显著性水平取 =0.05.H0: 2=02=1.5, H1: 2 02代码：X=49.4 50.5 50.7 51.7 49.8 47.9 49.2 51.4 48.9;var(X)h,p,muci,stats=vartest(X,1.5,0.05,both)运行结果：varX = 1.5278h = 0p = 0.8383muci = 0.6970 5.6072stats = chisqstat: 8.1481df: 8经检验，p 值=0.8383=0.05, 接受原假设
36、，故结论是：在显著水平=0.05 下，每包化肥质量的方差等于 1.5.4. 总体均值未知，两正态总体方差比较的 F 检验F 检验，又称方差齐性检验（方差齐性，即两方差没有明显差别） ，主要用于：均数差异的明显性检验、分别各有关因素并估量其对总变异的作用、剖析因素间的交互作用、方差齐性检验等。这里用来做方差齐性检验，代码：h,p,muci,stats=vartest2(X,Y, alpha, tail)例 16 还是例 14 中的样本数据，检验甲、乙两台机床加工的产品直径的方差是否相等？显著性水平取 =0.05.H0: 12=22, H1: 12 22代码：X=20.1 20.0 19.3 20.6 20.2 19.9 20.0 19.9 19.1 19.9;Y=18.6 19.1 20.0 20.0 20.0 19.7 19.9 19.6 20.2;h,p,muci,stats=vartest2(X,Y,0.05,both)运行结果：h = 0p = 0.5798muci = 0.1567 2.8001stats = fstat: 0.6826df1: 9df2: 8经检验，p 值=0.5798=0.05, 接受原假设，故结论是：在显著水平=0.05 下，两台机床加工的产品直径的方差没有显著差异。
