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逆函数的解析表达式,有些时候要做到这点比较困难,这限制了逆变换法的适用范围。
当无法给出分布函数F逆函数的解析表达式时,舍选法是另外的选择。舍选法的适用范围比逆变换法要大,只要给出概率密度函数的解析表达式即可,而大多数常用分布的概率密度函数是可以查到的。
算法思想:
给定轮廓,并在轮廓范围内生成均匀分布序列;
算法实现:
设概率密度函数为$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 = lemda*exp(-lemda*x);
ymax = lemda;
ymin = 0;
Y = (ymax-ymin)*rand(1,Len)-ymin;
data = x(Y<=fx);
结果图:
舍选法需要对数据重复筛选,性能不如逆变换法,但其适应性更广,无法得到分布函数的逆函数时,舍选法不失为一个选择。
其实本质就是取分布以内的随机数,颠倒过来,给定分布,将分布无限细分,生成各区间对应数量的随机数,也可以实现近似。
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$进行从0到∞的积分,容易得出$\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 = 2*pi*U2;
X = R.*cos(theta);
Y = R.*sin(theta);
效果图:
啰嗦一句:从scatter(X,Y)的结果中,可以直观看出X、Y不相关,但不能确定不独立。$X*Y = (X,Y)$与$P_X*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\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 = sigma*sqrt(-2*log(U1));
结果图:
例三:泊松分布(Poisson distribution)
背景介绍
一段时间内,事件发生的次数,组成计数过程。泊松分布是一种常用的计数过程。
生活中,大量事件是有固定频率的,即可以通过一段时间内计数找出规律:
某医院平均每小时出生3个婴儿
某公司平均每10分钟接到1个电话
某超市平均每天销售4包xx牌奶粉
某网站平均每分钟有2次访问
这里不展开讲泊松分布的来龙去脉,直接给出公式:
再来回顾指数分布,指数分布描述的是:计数过程中,相邻事件发生的时间间隔概率分布。
对应上面的问题,则有
婴儿出生的时间间隔
来电的时间间隔
奶粉销售的时间间隔
网站访问的时间间隔
指数分布可由泊松分布推出:
泊松分布 —> 指数分布:
设相邻两次事件间隔为$T$,起始时刻为$T_{start}$,则终止时刻为$T_{start}+T$,$P\{ T \ge t \}$表示$[T_{start},T_{start}+t]$时间内没有事件发生,即:
从而事件发生的概率为:
对应概率密度为:
对于泊松分布,得出结论:
两次事件的时间间隔为负指数分布,且均值为$\frac{1}{\lambda }$;
事件1与事件2,事件2与事件3,...事件n与事件n+1,相邻事件的时间间隔具有独立同分布的特性。
算法实现
此处给出的泊松分布随机数,产生原理基于指数分布,更多方法可以戳这里。
首先给出算法思想与算法步骤:
算法思想:
根据上文分析可知,Poisson分布对应一段时间内(记为$t_{max}$)事件发生次数,而指数分布exp对应相邻时间的时间分布;
反过来,已知两两相邻的时间变量,该变量服从指数分布且相互独立,所有时间变量相加—>并让其不超过一段时间的总量$t_{max}$,则累加数的分布对应泊松(Poisson)分布。
算法步骤:
生成一组均匀分布随机数:$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的结果对比):
二、离散型随机数
假设已经拥有U(0,1)的均匀分布数据,离散型随机数根据[0,1]划分区间,给出对应变量值即可。
例:离散随机变量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;
直方图结果:
可以看出结果服从指定的分布率。
补充:本文仅简要介绍一维分布随机数的生成方式,有两个主要的点没有涉及:
均匀分布随机数的产生方式;
多维随机数的产生方式。
参考:
matlab逆变换法产生随机数_信号处理——生成给定分布随机数相关推荐
- matlab生成不重复的随机数_怎么生成不重复随机数——《超级处理器》应用
生成随机数,大部分同学都会. 那么,如何生成,不重复的随机数呢?例如,怎么生成20个,100以内的不重复随机数?五秒时间,思考下怎么做? 问题挺简单,做起来还是比较复杂.如果用超级处理器,就非常方便, ...
- java生成指数分布随机数_生成特定分布随机数的方法
生成随机数是程序设计里常见的需求.一般的编程语言都会自带一个随机数生成函数,用于生成服从均匀分布的随机数.不过有时需要生成服从其它分布的随机数,例如高斯分布或指数分布等.有些编程语言已经有比较完善的实 ...
- 生成特定分布随机数的方法:Python seed() 函数numpy scikit-learn随机数据生成
描述 seed() 方法改变随机数生成器的种子,可以在调用其他随机模块函数之前调用此函数.. 语法 以下是 seed() 方法的语法: import random random.seed ( [x] ...
- 生成特定分布随机数的方法
生成随机数是程序设计里常见的需求.一般的编程语言都会自带一个随机数生成函数,用于生成服从均匀分布的随机数.不过有时需要生成服从其它分布的随机数,例如高斯分布或指数分布等.有些编程语言已经有比较完善的实 ...
- matlab逆变换法产生随机数_matlab数值积分方法(一)
一元函数定积分的数学表示为 在被积函数f(x)理论上不可积时,即无法求出该积分的解析解,所以往往要采用数值方法来求解.求解定积分的数值方法是多种多样的,如简单的梯形法.Simpson 法.Romber ...
- matlab逆变换法产生随机数_matlab中产生随机数的程序
1. 由 U ( 0,1 )分布的随机数产生 U ( a,b )的随机数 r=rand(1,20); s=a+(b-a)*r; 例: r=rand(1,20); s=2+(10-2)*r s = Co ...
- java 8位随机数_JAVA中生成指定位数随机数的方法总结
JAVA中生成指定位数随机数的方法很多,下面列举几种比较常用的方法. 方法一.通过Math类 1 public static String getRandom1(intlen) {2 int rs = ...
- 用一个随机数函数去生成另一个随机数函数:rand(a)生成rand(b)以及rand(a,b) 生成rand(c,d)
让我们把这个问题泛化一下,从特殊到一般.现在我给你两个生成随机数的函数Randa, Randb.Randa和Randb分别产生0到a的随机数和0到b的随机数,a,b不相等 (相等就没必要做转换了).现 ...
- php怎么不重复随机数,php怎么生成不重复随机数
php怎么生成不重复随机数 php生成不重复随机数的方法:首先利用range函数创建一个包含指定范围的元素的数组:然后利用shuffle函数把数组中的元素按随机顺序重新排列:最后取出数组中的一段元素即 ...
最新文章
- bzoj 1409 Password 矩阵快速幂+欧拉函数
- C#基础--类/接口/成员修饰符,多态、重载、重写,静态和非静态
- 手机下载Python_将安卓手机打造成 Python 全栈开发利器
- Python类私有方法的陷阱
- C语言 main函数
- 抖音做综艺,差点意思
- 了解Minimax算法
- 富龙热电:望眼欲穿矿难拿
- 嵌入式软件开发下的数据积累
- 了解arXiv,及arXiv的注册详细操作。
- 持续更新就是给软件上医保
- 方德系统装exe文件_国产处理器+自主OS完美运行exe程序?英特尔认为有侵权嫌疑...
- UML顺序图(sequence diagram)
- FL Studio教程之如何慢慢降音
- matlab巴特沃斯滤波器用法
- c语言实现lower_bound和upper_bound
- Python中的XML解析错误[Et.parse(xml) ‘gbk‘ codec can‘t decode byte]分析与解决
- 微信iPad协议(8013)
- 香港十大正规外汇黄金交易平台排名(2021版)
- python爬去网页数据并对比excel中数据是否一致_python入门之对比两份excel表格数据...
热门文章
- 模仿Airbnb的悬浮搜索框动画
- Effective Modern C++[实践]->优先使用nullptr,而非0或NULL
- 示例库 - 超过50个流程图 (Collection: Over 50 Flowchart Examples)
- 黄仁勋:串行计算过时并行计算是未来
- 刚子扯谈:未完待续的微信5.0
- 光学时钟“升天”助力卫星精准导航
- python+itk+读取dicom数据,并保存为nii文件
- 4年的数学竞赛学习,终于成就这位牛娃的北大保送之路!
- 最全的国外机器学习资源(上)
- Oracle - DBMS_LOB函数和用法