Matlab笔记——数值计算—概率篇017

17. 数值计算—概率篇

一、计算组合数、排列数

!n——factorial(n)或prod(1:n)

k

C——nchoosek(n,k)

n

k

A——factorial(n)/factorial(n-k)

n

二、生成随机数

1. rand(m,n)

——生成m×n的服从[0,1]上均匀分布的随机数；

用a + (b-a).*rand(m,n)生成m×n的服从[a,b]上均匀分布的随机数。

2. 二项分布与正态分布随机数

binornd(N,P,m,n)——生成m×n的服从二项分布B(N,P)的随机数；

normrnd(MU,SIGMA,m,n)

——生成m×n的服从正态分布N(MU,SIGMA2)的随机数；

3. 通用格式：

分布缩写+rnd(分布参数, m,n)

或random(‘分布名或缩写’, 分布参数, m,n)

可以用来生成m×n该分布的随机数。各种分布名见下图：

4. 使用randsample和randsrc函数生成指定离散分布随机数

X=randsample(N, k, replace, w)

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

X= randsrc(m,n,[x; p])

生成m×n的随机矩阵，服从取值为向量x, 对应概率为向量p的离散分布。

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

生成服从3×5的该分布的随机数。

代码：

xvalue = [-2 -1 0 1 2];

xp = [0.05 0.2 0.5 0.2 0.05];

% 调用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 0

0 0 0 -1 -1

1 1 0 0 1

y = -1 -1 1 1 -1

-1 0 0 2 0

-1 0 -1 0 0

5. 已知概率密度函数，生成服从该分布的随机数 例2 设随机变量X 的概率密度函数为（抛物线分布）：

6(1), 01

() 0, x x x f x -<

调用crnd 函数（来自《MATLAB 统计分析与应用 40个案例分析》作者：谢中华），生成3×5个服从该分布的随机数。

代码：

pdffun = '6*x*(1-x)'; % 密度函数表达式 x = crnd(pdffun,[0 1],3,5)

运行结果：x = 0.3160 0.6866 0.2724 0.2816 0.1268 0.2681 0.8439 0.1948 0.7999 0.5383 0.7377 0.2040 0.4932 0.1948 0.6909

6. 生成多元分布的随机数

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组的二元正态分布随机数，其中分布的参数为

1013==20316μ????

∑ ? ?????

， 代码：

mu = [10 20];

sigma = [1 3; 3 16];

xy = mvnrnd(mu, sigma, 3)

运行结果：xy = 11.8336 25.7385

9.0347 17.8026

9.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绘制卡方分布密度函数在自由度分别为1、5、15的图形。

代码：

x=0:0.1:30;

y1=chi2pdf(x,1); plot(x,y1,':')

hold on

y2=chi2pdf(x,3);plot(x,y2,'+')

y3=chi2pdf(x,10);plot(x,y3,'o')

axis([0,30,0,0.2])

运行结果：

四、随机变量的分布函数

分布函数定义为：F(x)=P{X≤x}，表示随机变量X的取值落在(-∞,x)范围内的概率。引入分布函数的目的，就是可以计算随机变量X的取值落在任意区间内的概率，例如，

P{a

1. 通用格式：

cdf(‘分布名或缩写’, x, 分布参数)

——返回该分布的分布函数；

例如，Pk=cdf('bino',3, 10, 0.4)

2. 专用函数

分布名缩写+cdf(x, 分布参数)

例如，binocdf (k, n, p)

五、逆分布函数

已知F(x)=P{X≤x}的值，求x点。

1. 通用格式：

icdf(‘分布名或缩写’, p, 分布参数)

——返回该分布的分布函数F(x)=P{X≤x}=p的x值；

例如，Pk=icdf('bino',3, 10, 0.4)

2. 专用函数

分布名缩写+inv(p, 分布参数)

例如，binoinv (p, N, P)

3. 分位数

N(0,1)的上α分位数：P{X>uα}=α，故

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) 几何平均值

11n

n

i i G x =??=∏ ???

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

代码：geomean(X) (4) 调和平均值

11

n

i i

n H x

==

∑

即先取倒数，再取算术平均值，再取倒数回来。例如前半段时速60公里，后半段时速30公里（两段距离相等），则其平均速度为两者的调和平均数时速40公里。在实际中，往往由于缺乏总体单位数的资料而不能直接计算算术平均数，这时需用调和平均法来求得平均数。

代码：harmmean(X) (5) 众数

出现次数最多的数，代码：mode(X)

2. 期望和方差 (1)样本数据

mean(X)——样本均值1

1n

i i x x n ==∑

var(X)或var(X,0)——样本方差2

21

1=()1n

i i S x x n =--∑ var(X,1)——方差21

1()n

i i D x x n ==-∑

std(X)或std(X,0)——样本标准差S std(X,1)

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

[E,D]=分布名缩写+stat(分布参数)

返回E 为期望，D 为方差；例如，[E,D]=normstat(MU,SIGMA)

(3) 已知离散型随机变量的分布律，求期望和方差 例5设离散型随机变量X 服从如下分布：

求E(X), E(X 2-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*p’ Y=X.^2-1;

EY=sum(Y.*p)

XX=X.^2;

EXX=sum(XX.*p);

DX=EXX-EX^2

运行结果：EX = 0EY = 1.6000DX = 2.6000

3. 极差、偏度、峰度

极差，数据的最大值与最小值之差。代码：range(X)

偏度，是数据关于均值不对称的指标。若偏度为负，说明均值左边的数据比右边的数据更散；若偏度为正，说明均值右边的数据比左边的数据更散，正态分布的偏度为0. 代码：skewness(X) 峰度，是用来反映数据的分布曲线顶端尖峭或扁平程度的指标。代码：kurtosis(X)

4. 矩

k阶中心矩——m=moment(X, k)

k阶原点矩——可以使用下面自编的OriginMoment.m

function y=OriginMoment(X,k);

% X为样本矩阵

[n,m]=size(X);

y=zeros(1,m);

for ii=1:m;

y(ii)=sum(X(:,ii).^k)/n;

end

5. 协方差（矩阵）、相关系数（矩阵）

对于样本数据矩阵：

[]1112121

2221

21

2

n n n m m mn m n

x x x x x x X X X X x x x ?????

??==??????

Matlab 将矩阵X 的每一列X j , j=1,…,n 作为一个随机变量的样本，每一行[x j1,x j2, …, x jn ], j=1,2,...,m 作为n 个随机变量的联合分布的一个样本。

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

样本协方差公式：

1

1(,)()()1n

i i i COV X Y x x y y n ==---∑ 111212122212n n n n nn n n

c c c c c c c

c c ?????

??∑=

?????? , 其中(,)ij i j

c COV X X =

称为随机变量X 1, …, X n 的协方差矩阵（实对称、非负定）。

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

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

代码：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.8444 9.1065

0.0463 3.9978 1.8185

7.7491 2.5987 2.6380

8.1730 8.0007 1.4554

CovX = 19.6061 -6.3812 5.3900

-6.3812 11.6214 -5.5894

5.3900 -5.5894 9.7937

Cov12 = 19.6061 -6.3812

-6.3812 11.6214

Cov11 = 19.6061 19.6061

19.6061 19.6061

VarX1 = 19.6061

RhoX = 1.0000 -0.4227 0.3890

-0.4227 1.0000 -0.5239 0.3890 -0.5239 1.0000

七、频率直方图

1. 统计数值或字符出现的频数、频率

设X 为数值型或字符型向量，代码：tabulate(X); 输出X 的值、频数、频率的三列表。

注意：若X 中的值都是非负整数，则min(X)与max(X)之间的任一整数都统计，未出现的频数=0.

2. 经验分布函数

1,, ()n n x x x F x n

≤=

中的个数

实际上是累积频率直方图曲线。依据样本以频率估计概率，当n 充分大时，()n F x 是总体分布函数()F x 的近似。

例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') % 绘制经验分布函数图形(连线型） hold on

[h,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.9801

mean: 3.2968 median: 2.1667 std: 3.2969

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; % 分组宽度

x

F (x )

Empirical CDF

[r,xout]=hist(X,[a:d:b]); % 或者[r,xout]=hist(X,m);

% 计算直方图数据, r返回每段的频数, xout返回每段的中心位置

f=r./(n*d); % 计算频率

bar(xout,f); % 绘制频率直方图

hold on

x=0:0.01:10;

y=fpdf(x,3,5);

plot(x,y,'k-');

axis([0 10 0 1]);

title('频率密度直方图');

运行结果：

频率密度直方图

0246810注：也可以用例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)

运行结果：

八、箱线图

又称为盒形图。在一条数轴上，以数据的上下四分位数（Q 1-Q 3）为界画一个矩形盒子（中间50%的数据落在盒内）；在数据的中位数

1

2

3

4

56

7

8

9

0.003

0.010.020.050.100.250.500.75

0.900.950.980.99

0.997Data

P r o b a b i l i t y

Normal Probability Plot

位置画一条线段为中位线；默认延长线不超过盒长的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);

x = [x1 x2];

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)

对于二项分布为：[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.08

15.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.1397

muci = 14.9561

15.1559

sigmaci = 0.0961

0.2550

logL = -5.9933

info = 0.0020 0.0000

0.0000 0.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值在不在阴影部分（拒绝域），空白部分为接受域。