1.9 概率统计
概率统计属于数学上的一个大的分支，用于统计自然界随机现象发生的规律，概率统计应用较广泛，例如股票的预测、列车客流及其他信息统计等。本节着重讲述概率密度函数、随机变量特征分析、随机数概率密度函数绘图以及蒙特卡诺方法等。
1.9.1 概率密度函数
MATLAB统计工具箱是一套完备的统计分析软件，表1-13给出了MATLAB概率分布工具箱函数，读者可以根据表1-13所示的分布函数，借助MATLAB函数帮助功能进行相应问题的求解。
表1-13 MATLAB概率分布工具箱函数
通过相应的MATLAB中的求概率密度函数，也可以求表1-13中所示的23种输入分布的概率密度函数，此函数的调用格式如下：
Y Y pdf('name', X, A1, A2, A3)
其中，name为表1-13中分布的字母缩写，X为样本矩阵，A1、A2和A3是分布参数矩阵，Y为概率密度矩阵。对于某些分布，有些参数可以不必输入。X、A1、A2和A3必须是具有同样大小的矩阵，当其中之一输入为标量时，程序会将其调整为与其他输入同维的矩阵。
具体如表1-14所示。
表1-14 概率密度函数使用格式
1.9.2 随机变量特征分析
1．期望
设离散性随机变量的分布律为：
则相应X的期望如下：
而对于来自总体X的一个样本，设其样本值为
，则定义样本均值为：
则
依概率收敛于X的均值。
在MATLAB统计工具箱中，提供了求解随机变量均值的函数mean()，具体如下：
>> mean(rand(5,1))
ans =
0.5440
2．方差、标准差、矩
方差是用来刻画随机变量X取值分散程度的一个量，其一般用下式表达：
在应用上还引入与随机变量X具有相同量纲的量
，记为σ（X），称为标准差或均方差。
X的k阶中心矩应为：
可知方差即为二阶中心矩。
对于一个样本来说，样本方差通常分为无偏估计和有偏估计，具体如下。
无偏估计式：
有偏估计式：
样本标准差也对应有如下两种形式：
或
样本的k阶中心矩为：
MATLAB工具箱中提供了可供用户求解的方差函数var()、标准差函数std()和矩函数moment()，具体如下：
>> var(rand(5,1))
ans =
0.0799
>> std(rand(5,1))
ans =
0.3791
>> moment(rand(5,1),1)
ans =
0
3．协方差、相关系数
随机变量X、Y的协方差和相关系数的定义式如下：
对于n维随机变量，通常用协方差矩阵描述它的二阶中心矩。如对于二维随机变量(X,Y)，定义协方差矩阵形式为：
其中：
MATLAB提供了求解协方差和相关系数的函数：cov()函数和corrcoef()函数，具体如下：
>> cov(rand(5,2))
ans =
0.1607 -0.0266
-0.0266 0.0536
>> corrcoef(rand(5,2))
ans =
1.0000 0.1471
0.1471 1.0000
1.9.3 随机数概率密度函数绘图
MATLAB统计学工具箱，自带很多产生服从某个函数分布的随机数的函数，用户可以很方便地调用和使用，具体的介绍如下。
（1）均匀分布
生成(0,1)区间上均匀分布的随机变量，MATLAB编程如下：
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
x=rand(100000,1); % 100000行1列
hist(x,50);
运行程序可得到图形如图1-44所示。
图1-44 均匀分布
（2）正态分布
生成服从标准正态分布（均值为0，方差为1）的随机数，MATLAB编程如下：
x=randn(10000,1); % 100000行1列
hist(x,50);
运行程序可得到图形如图1-45所示。
图1-45 正态分布
（3）卡方（Chi-square）分布
卡方分布只有一个参数：自由度v，MATLAB编程如下：
x=chi2rnd(5,100000,1); %卡方分布的自由度是5，100000行1列
hist(x,50);
运行程序可得到图形如图1-46所示。
图1-46 卡方分布
（4）F分布
F分布有两个参数：v1和v2，MATLAB编程如下：
x=frnd(3,5,1000,1); % (v1=3,v2=5)的F分布
hist(x,50);
运行程序可得图形如图1-47所示。
图1-47 F分布
（5）t分布
t分布有1个参数：自由度v。MATLAB函数调用如下：
x=trnd(7,100000,1); % 参数为(v=7)的t分布
hist(x,50);
运行程序可得到图形如图1-48所示。
图1-48 t分布
（6）Beta分布
Beta分布有两个参数，分别是A和B，MATLAB函数调用如下：
x=betarnd(2,5,100000,1); % 产生A=2，B=5的Beta分布
hist(x,50);
运行程序可得到图形如图1-49所示。
图1-49 beta分布
（7）指数分布
指数分布只有一个参数：mu，生成指数分布随机数的MATLAB函数调用如下：
x=exprnd(3,100000,1); % 产生mu=3的指数分布
hist(x,50);
运行程序可得到图形如图1-50所示。
图1-50 指数分布
（8）Gamma分布
Gamma分布有两个参数：A和B，生成Gamma分布随机数的MATLAB函数调用如下：
x=gamrnd(2,5,100000,1); % 产生A=2，B=5的Gamma分布
hist(x,50);
运行程序可得到图形如图1-51所示。
图1-51 Gamma分布
（9）对数正态分布
对数正态分布有两个参数：mu和sigma，生成对数正态分布随机数的MATLAB函数调用如下：
x=lognrnd(-1,0.5,1000,1); % 产生mu=-1，sigma=0.5的对数正态分布
hist(x,50);
运行程序可得到图形如图1-52所示。
图1-52 对数正态分布
（10）瑞利（Rayleigh）分布
瑞利（Rayleigh）分布有1个参数：B。生成瑞利分布随机数的MATLAB函数调用如下：
x=lognrnd(-1,0.5,1000,1); % 产生B=2的瑞利（Rayleigh）分布
hist(x,50); % 直方图
运行程序可得到图形如图1-53所示。
图1-53 瑞利分布
（11）威布尔（Weibull）分布
威布尔（Weibull）分布有两个参数：scale参数A和shape参数B。生成Weibull分布随机数的MATLAB函数调用如下：
x=lognrnd(-1,0.5,1000,1); % 产生A=3，B=2的威布尔分布
hist(x,50); % 直方图
运行程序可得到图形如图1-54所示。
图1-54 威布尔分布
（12）二项分布
二项分布有两个参数：n和p。MATLAB编程如下：
x=binornd(10,0.45,100000,1); % n=10，p=0.45
hist(x,11); % 直方图
运行程序可得到图形如图1-55所示。
图1-55 二项分布
（13）几何分布
几何分布的参数只有1个：p，MATLAB编程如下：
x=geornd(0.4,100000,1); % p=0.4
hist(x,50); % 直方图
运行程序可得到图形如图1-56所示。
图1-56 几何分布
（14）泊松（Poisson）分布
泊松（Poisson）分布的参数只有1个：lambda，此参数要大于零。MATLAB函数调用如下：
x=poissrnd(10,100000,1); % lambda = 10
hist(x,50); % 直方图
运行程序可得到图形如图1-57所示。
图1-57 泊松分布
1.9.4 蒙特卡洛Monte Carlo算法
Monte Carlo算法计算的结果收敛的理论依据来自于大数定律，Monte Carlo算法计算要进行很多次抽样，才会比较好地显示出来。如果Monte Carlo计算结果的某些高阶距存在，即使抽样数量不太多，这些渐进属性也可以很快地达到。
在金融产品定价中，大多数问题求解基于某个随机变量的函数的期望值。考虑一个欧式期权，假定已知在期权行权日的股票服从某种分布（理论模型中一般是正态分布），那么用期权收益在这种分布上做积分求期望值即可。
给定曲线
和曲线
，曲线的交点为：
、
。曲线围成平面有限区域，用蒙特卡洛方法计算区域面积，编程如下：
% Designed by Yu Shengwei From SWJTU University
% 2014年12月29日
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
format short
tic % 运算计时
P=rand(10000,2);
x=2*P(:,1)-1;
y=2*P(:,2);
points=find(y<=2-x.^2&y.^3>=x.^2);
M=length(points);
S=4*M/10000
figure('color',[1,1,1])
plot(x(points),y(points),'bs')
toc % 计时结束
运行程序输出结果如下：
S =
2.1168
时间已过 0.085066 秒。
得到相应的图形如图1-58所示。
图1-58 面积求解
蒙特卡洛方法计算冰淇淋锥体积
，
，又冰淇淋锥含于体积为8的六面体
；N个点均匀分布于六面体中，锥体中占有m个，则锥体与六面体体积之比近似为m∶N，即
。
绘制该冰淇淋锥体，编程如下：
% Designed by Yu Shengwei From SWJTU University
% 2014年12月29日
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
format short
tic % 运算计时
m=20; % 六面体中点数
n=100; % 锥体中点数
t=linspace(0,2*pi,n); % 点
r=linspace(0,1,m); % 半径
x=r'*cos(t); % x
y=r'*sin(t); % y
z1=sqrt(x.^2+y.^2); % 半径
z2=1+sqrt(1+eps-x.^2-y.^2); % 体积2
X=[x;x];Y=[y;y]; % 矩阵扩充
Z=[z1;z2]; % 矩阵扩充
figure('color',[1,1,1])
mesh(X,Y,Z) % 3D曲面
view(0,-18) % 视角
colormap([0 0 1]),axis off % 取消轴显示
toc % 计时结束
运行程序输出结果如图1-59所示。
图1-59 冰淇淋模型
如图1-59所示冰淇淋模型，该冰淇淋模型体积求解如下：
设
，
，
，则该体积转化为一个定积分表达式如下：
MATLAB编程如下：
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
format short
tic % 运算计时
L=7; % 重复计算次数
N=10000; % 点数
for k=1:L
P=rand(N,3); % 随机数
x=2*P(:,1)-1; % x
y=2*P(:,2)-1; % y
z=2*P(:,3); % z
R2=x.^2+y.^2; % 半径平方
R=sqrt(R2); % 半径
II=find(z>=R&z<=1+sqrt(1-R2)); % 满足不等式
m=length(II);
q(k)=8*m/N; % 体积
end
error=q-pi % 误差
figure('color',[1,1,1])
plot(x(II),y(II),'bp') % 画出x,y视图
toc % 计时结束
运行程序输出结果如下：
error =
Columns 1 through 5
0.0096 -0.0336 -0.0600 0.0008 0.0488
Columns 6 through 7
-0.0344 0.0080
输出图形如图1-60所示。
图1-60 xy平面映射