matlab逆变换法产生随机数_信号处理——生成给定分布随机数

作者：桂。

时间：2017-03-12 19:31:55

前言

该文作为曲线拟合与分布拟合一文的引言部分，在进行分布拟合的过程中遇到了一个问题：如何确定分布的随机数？该文主要讨论了这一问题。

1)连续型随机数；

2)离散型随机数；

本文内容为自己的学习笔记，内容多有借鉴他人，在最后一并给出链接。

一、连续型随机数

假设已经拥有U(0，1)的均匀分布数据。

A-逆变换法(Inverse Transform Method)

若我们能提供概率分布F及其逆函数的显式形式，则可通过简便方法生成所需随机数。

性质：

U是服从[0，1]均匀分布的随机变量，如果X = F^{-1}(U)，则X的分布函数与F相同，即F_X(x) = F.

证明：

F_X(a) = P(X \le a) = P(F^{-1}(U) \le a) = P(U \le F(a)) = F(a).即F_X与F相同，得证。

算法步骤：

生成一个服从均匀分布的随机数U~Unit(0,1);

设F(x)为指定分布的分布函数，则X = F^{-1}(U)即为指定分布的随机数。

示例：生成满足\lambda = 2的指数分布随机数。

分析：

由f(x)得出F(x) —>F(x) = 1 - {e^{ - \lambda x}}，进而求得F(x)逆函数，得出X = {F^{ - 1}}(u) = - \frac{1}{\lambda }\ln (1 - u).

代码：

Len = 1000000;

u = rand(1,Len);

lemda = 2;

x = -1/lemda*(log(1-u));

对应结果：

B-舍选法(Acceptance-Rejection Method)

通常情况下

在无法直接获得F^{-1}(x)的概率分布函数解析表达式的情况下，在实际应用中另一种可选的方法是舍选法。与逆变换方法相比，在满足一定条件下该方法的应用范围更为广泛——即仅需提供相应的概率密度函数解析表达式即可实现模拟抽样操作；而这种抽样技术尤其适用于那些概率质量函数或概率密度函数具有明确形式但累积分布函数却难以求解的情形——例如许多常用分布在统计应用中较为常见但其累积分布函数却不易直接计算或不存在闭合形式解的情况。

算法思想：

给定轮廓，并在轮廓范围内生成均匀分布序列；

算法实现：

设概率密度函数为f(x)，首先生成一个均匀分布随机数X~Unit(x_{min}，x_{max});

独立地生成一个均匀分布随机数Y~Unit(y_{min}，y_{max});

若Y \le f(X)，则返回X，否则重复第一步。

给出一张舍选法的原理图：

代码(data即为最终的随机数)：

xmin = 0;

xmax = 5;

Len = 1000000;

x = (xmax-xmin)*rand(1,Len)-xmin;

lemda = 2;

fx = lemdaexp(-lemdax);

ymax = lemda;

ymin = 0;

Y = (ymax-ymin)*rand(1,Len)-ymin;

data = x(Y<=fx);

结果图：

舍选法在应用时需要反复进行筛选操作，在性能方面相对而言 inverse transform method 具有更高的效率。然而，在适应性方面却更为广泛。当无法获得该分布函数的反函数时，在某些情况下仍可作为一种可行的方法。

本质上就是抽取分布内部的随机数，并将其反转过来。给定一个特定的分布函数，在不对其概率密度函数求导的情况下，将该分布函数无限细分区间，并在每个区间内生成相应数量的随机数也是一种近似方法。
同样也能实现近似。

C-组合法

当目标分布可以通过其他分布经过基本算术运算表示时, 可以使用组合算法来产生相应的随机数. 本节将通过几个具体的例子来简要阐述这一过程.

例一：正态分布(Box Muller方法)

生成正态分布随机数的方法不止一种，在此之外的经典方法是Box-Muller方法。具体而言，任何正态分布都可通过相应的标准正态分布进行转换以实现生成。

Box Muller理论推导：

令(X,Y)为两个互相独立且服从标准正态分布N(0,1)的随机变量，则其概率密度函数具有指数形式缺少一个负号。

转化为极坐标形式，令

，则R有分布：

Eq.(1)

令

，则分布函数的反函数为：

由逆变换法性质可知：

满足F_R的分布可由1-Z_{$U(0,1)$得出，亦即：可由$Z$}U(0,1)得出；

进一步分析P_{\theta}的方法，则同样地，在Eq.(1)中对该表达式关于r从零到无穷大的积分进行计算后，则可推导出\theta服从均匀分布；由此可知，该随机数可以直接通过均匀分布生成。

显然地，我们有P(r;\theta)=P(r)P(\theta)成立；这表明二者相互独立；由此可知, 它们可以通过各自的均匀分布独立生成.

以上为Box Muller的原理分析。

Box Muller算法实现：

分别生成两组均匀分布随机数：U_1_{$U(0,1)$、$U_2$}U(0,1)；

生成R以及\theta的分布：

生成两组独立的标准正态分布:

生成符合给定均值、方差具体要求的正态分布；

代码(生成标准正态分布)：

Len = 10000;

U1 = rand(1,Len);

U2 = rand(1,Len);

R = sqrt(-2*log(U1));

theta = 2piU2;

X = R.*cos(theta);

Y = R.*sin(theta);

效果图：

简而言之：在scatter(X, Y)的结果中可以看出X、Y线性不相关但无法断定它们之间不存在非线性关系.x \cdot y = (x,y)与P_x \cdot P_y = P(x,y)不是一个概念。同一个圆形区域中X、Y可以以不同的密度分布存在于其中。

再看看均匀分布的scatter(X,Y)：

scatter结果：一个圆形、一个方形，但两种情况下X、Y都是独立同分布的。

例二：瑞利分布(Rayleigh Distribution)

其实Box Muller方法的中间过程，包含了瑞利分布，即P(r)。给出PDF定义式：

f(r) = \frac{r}{{{\sigma ^2}}}{e^{ - \frac{{{r^2}}}{{2{\sigma ^2}}}}}

其中\sigma > 0.

理论分析：

令(X,Y)为一对相互独立地遵循正态分布N(0,\sigma ^2)的随机变量，则其概率密度函数为f(x,y)。

f\left( {x,y} \right) = f(x)f(y) = \frac{1}{{2\pi {\sigma ^2}}}{e^{ - \frac{{{x^2} + {y^2}}}{{2{\sigma ^2}}}}}

转化为极坐标，dxdy = rdrd\theta ，并对\theta求取积分，得：

f\left( r \right) = \frac{r}{{{\sigma^2}}}{e^{ - \frac{{{r^2}}}{{2{\sigma ^2}}}}}

假设随机变量(X,Y)是相互独立且均服从标准正态分布N(0,\sigma^2)的一对变量，则通过Box-Muller方法生成瑞利变量R即可得到所要求的瑞利分布。

回顾公式：R = \sqrt {{X^2} + {Y^2}} ，两个独立同正态分布(零均值)信号的包络是瑞利分布。

算法步骤：

生成一组均匀分布随机数：U~U(0,1)；

根据给定的参数\sigma，生成：R = \sigma \sqrt { - 2\ln U} ;R即为满足要求的瑞利分布；

对应代码：

Len = 1000000;

U1 = rand(1,Len);

U2 = rand(1,Len);

sigma = 2;

R = sigmasqrt(-2log(U1));

结果图：

例三：泊松分布(Poisson distribution)

背景介绍

在某个时间段内，发生次数构成了计数过程。泊松分布是一种典型的计数过程。

生活中，大量事件是有固定频率的，即可以通过一段时间内计数找出规律：

某医院平均每小时出生3个婴儿

某公司平均每10分钟接到1个电话

某超市平均每天销售4包xx牌奶粉

某网站平均每分钟有2次访问

这里不展开讲泊松分布的来龙去脉，直接给出公式：

再次学习关于指数分布的知识：在计数过程中的相邻事件发生的时间间隔的概率情况。

对应上面的问题，则有

婴儿出生的时间间隔

来电的时间间隔

奶粉销售的时间间隔

网站访问的时间间隔

指数分布可由泊松分布推出：

泊松分布 —> 指数分布：

设两次相邻事件的时间间隔为随机变量T，将其初始时刻定义为t_0=T_{start}。那么结束时刻则表示为初始时刻加上时间间隔t_0+T=t_start+ T= T_start + T. 概率\text{P}\{ T \ge t \}代表在时间段\left[ t_0, t_0 + t \right]内没有任何事件发生。

从而事件发生的概率为：

对应概率密度为：

对于泊松分布，得出结论：

两次事件的时间间隔为负指数分布，且均值为\frac{1}{\lambda };

每个连续的两个相邻事件之间存在某种关系，并且这些时间间隔都满足独立同分布的特性

算法实现

我们在此处提供的服从泊松分布的随机数样本，在其生成机制建立在指数分布的基础上进行计算与模拟。如果您想了解更多的生成方法，请点击原文链接进一步了解相关技术细节。

首先给出算法思想与算法步骤：

算法思想：

基于上文的分析结果表明，在时间段内（以t_{max}计），Poisson分布代表事件的发生次数；而指数分布在描述相邻时间间隔方面起着关键作用。

反过来思考，在两个相邻时间段内考虑随机发生的事件次数，并假设这些时间间隔都服从指数分布且相互独立。进一步地，并要求其总和不超过某个最大时间间隔T_{max}。那么，在这种情况下累积事件数目遵循泊松过程的概率规律。

算法步骤：

生成一组均匀分布随机数：U~U(0,1);

利用逆变换法生成一系列独立的指数分布X_i ~ exp(\lambda );

记

如果Y > t_{max}，则停止，并输出k-1;否则，继续生成X_{k+1}，直到Y > t_{max}为止；

循环操作过程3；

输出的一系列整数(值为k-1)服从参数为\mu = \lambda t_{max}的Poisson分布。

考虑到上文已经详细说明了如何生成指数分布，在代码中我们直接调用指数随机量exprnd函数，并无需额外生成；本文着重讨论思路，在此过程中我们列举的例子均基于现成的函数调用

给出代码(data即为生成的数据)：

Iter = 10000; % 迭代次数，即data个数

lambda = 2;

t_max = 11;

%phei = lambda*t_max;%均值

for m=1:Iter,

k=1;

Y(m,1)=exprnd(1/lambda);

while Y(m,k) <= t_max % 时间间隔 (0,t_max]

Y(m,k+1)=Y(m,k)+exprnd(1/lambda);

k=k+1;

end

data(m)=k-1;

end

测试结果(将代码生成结果与MATLAB自带函数poissrnd.m的结果对比)：

二、离散型随机数

假设具有均匀分布在区间内的数据。将该区间划分为若干子区间后，在每个子区间内生成对应的随机数值。

例：离散随机变量X的分布率如下表，试生成该随机数。

X

-1

1

3

P(x)

0.7

0.2

0.1

tag = rand(1,1000);

data = tag;

data(tag>0.9) = 3;

data(tag<=0.7) = -1;

data(abs(data)<1) = 1;

直方图结果：

可以看出结果服从指定的分布率。

补充说明的是，在本文中我们仅仅进行了概述性介绍一维分布随机数的生成方式，并且未涵盖两个关键方面

均匀分布随机数的产生方式；

多维随机数的产生方式。

参考：