1.9 概率统计_MATLAB GUI设计入门与实战-QQ阅读中文仙侠网

书名：MATLAB GUI设计入门与实战
作者名：余胜威吴婷罗建桥编著
本章字数：2737字
更新时间：2021-03-12 17:50:22

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平面映射