基于Matlab语言的Monte Carlo入门教程(6)

（六）特殊连续分布

这里我将Matlab中没有对应函数的分布称为特殊分布。有多种方法可以用于生产服从这些分布的随机数。这里主要介绍两种最常见的。

1．逆CDF函数法

如果我们已知某特定一维分布的CDF函数，经过如下几个步骤即可生成符合该分布的随机数。（其中数学推导等在此处略去，详见相关数学书籍）

1. 计算CDF函数的反函数：

2. 生成服从(0,1)区间上均匀分布的初始随机数a

3. 令x=，则x即服从我们需要的特定分布的随机数。

为了更形象解说这种方法，这里选取柯西（Cauchy）分布作为例子。有时也称其为洛仑兹分布或者 Breit-Wigner 分布。柯西分布有一大特点就是，它是肥尾（Fat-tail，又译作胖尾）分布。在金融市场中，肥尾分布越来越受到重视，因为在传统的正态分布基本不考虑像当前次贷危机等极端情况，而肥尾分布则能很好地将很极端的情形考虑进去。 上图是Cauchy 分布和标准正态分布PDF图对比，看看是不是Cauchy分布的尾巴（x轴两端）更“胖”一点？

柯西分布的PDF函数是：

简化起见我们只考虑x0=0, γ=1情形。此时PDF函数是：

PDF函数对x作积分，就得到CDF函数（推导过程略）：

现在我们套用这三个步骤来生成服从Cauchy分布的随机数：

1. 计算得到Cauchy分布CDF函数的反函数为：

2. 使用rand()函数生成（0，1）区间上均匀分布的初始随机数。我习惯一次生成一堆这种随机数。

original_x=rand(1,100000);

3. 将初始随机数代入CDF反函数即可得到我们需要的Cauchy随机数。。 cauchy_x=tan((original_x-1/2)*pi);

上面这两句代码结合起来就生成了10万个服从参数为(x0=0, γ=1)Cauchy分布的随机数。

这种方法生成随机数与Cauchy分布有多大相似之处呢？这里有一个图可以说明： 蓝色的图形就是用hist画出的随机数的样本分布情况，红色线条是Cauchy分布理论的PDF图形。由此可看出生成的随机数挺符合Cauchy分布。

注意：上图中，我略去了x轴小于-12.5和大于12.5部分的图形——因为Cauchy是胖尾分布，会生成出的不少取值很大的随机数，而那些很大的值使得我们不方便用hist函数来画随机数分布图。

注意，这种方法本身虽然很简单，效率也很高，但有如下受限之处：

1．它有个可能会出错的地方，有的CDF函数的反函数在0或者1处的值是正/负无穷，例如此处的Cauchy分布就是这样，倘若用(0,1)均匀分布产生的初始随机数中包含0或者1，那么这个程序会出错。幸运的是，迄今为止，我用Matlab的rand()函数生成的随机数中还没有出现过0或者1。但不同版本的Matlab的这种情况也许会改变。此处提醒你，如果程序出错，不要忘记检查是否是这个错误。

2．CDF函数必须严格单调递增，这也就意味着，PDF函数在x定义域内必须处处严格大于零，否则CDF的反函数不存在。

3．即使CDF函数存在，如果它太复杂，可能导致计算速度太慢，甚至无法计算的后果。

2．接受/拒绝法 Acceptance-Rejection Method

Accelptence-Rejection方法的精髓在于“形似”，可以形象地将其比喻为制作冰雕——二者相同之处在于都要首先堆砌出雏形，然后再用将多出的部分削去。用此法生成服从f(x)分布的随机数，分为如下几大步骤：

1．首先，选用某个分布，如pdf为g(x)的分布，此时要计算一个常数c, 使得,对x定义域内任意的x都成立——这相当于使cg(x)图形完全“覆盖”住f(x)图形，容易理解，做冰雕时，最初堆出来的那堆冰块要比最终得到的雕塑大。

2．生成服从pdf为g(x)分布的随机数，假设生成的随机数为x0。 3．再生成一个服从（0，1）间的均匀分布的随机数y 4．如果，丢弃生成的x0；反之，生成的x0就是我们需要的、服从f(x)分布的随机数。

下面用一个例子结合图形解释这种方法，假设我们要生成的分布是:

，此pdf图形如下图的蓝色曲线。

1．我们选用（0，2）之间的均匀分布作为原始分布，即g(x)=0.5，此分布的pdf图见

下图中的绿色线。由条件：无论哪个x，都要成立，我们计算得到c要大于等于10.8。这种情况下，我们一般选择c=1.875。因为c选得越大，意味着我们堆砌的原始雏形越大，需要削去的部分越多，效率越低，所以我们要使得c尽量地小。

2．生成服从（0，2）之间的均匀分布的随机数，设它为x0

X0=unifrnd(0,2);

3．然后再生成一个服从（0，1）间的均匀分布的随机数y Y=rand;

4．如果，丢弃生成的x0，重新生成；反之，生成的x0就是我们需要的、服从f(x)分布的随机数，用于做后续计算。 以上步骤每次只能处理一个随机数，效率较低，下面这段代码可以一次性生成一堆随机数。

N=400000;c=1.875;gx=0.5 x0=unifrnd(0,2,1,N); y=rand(1,N);

fx0=(x0-0.5).*(x0-0.5)/2.4;

final_x=x0(y<=fx0./c/gx);

在视频教程中我会逐句解释每句含义，如果没听懂，一般是因为你对Matlab向量运算不熟悉，请参照Matlab基础教程学习此部分的内容，后面章节会有很多地方用得上。

这段语句生成的变量final_x即为服从f(x)分布的随机数组成的一个行向量。我们可以用hist查看这些随机数大致的分布。 hist(final_x,50);title('f(x)=(x-0.5)^2/2.4'); 如图所示，生成的随机数挺符合f(x)分布。

这种方法很简单，也不需要计算CDF函数的反函数，但它也有如下受限之处： 1．由于我们用随机数y来控制是否削去某个随机数x0，所以我们无法准确预知最终得到的随机数数量多少。

2．选择合适的g(x)分布是此方法最关键的技巧所在。g(x)的选择原则是在完全覆盖f(x)的前提下尽可能与f(x)形似，二者形状越相似，需要削去的部分就越少，这种方法的效率就越高。需要记住的是：很多时候，人们不选用这种方法的原因几乎都在于它的效率过低。

（七）特殊离散分布

离散分布关键在获得它的分布律，有了分布律我们计算骰子掷出点数小于等于某个数字的累积概率分布。一个简单的例子，假设我们有一个不均匀的骰子，获得六个点数的概率分别是：

点数 1 2 3 4 5 6 概率 .1 累积点数 1 累积概率 .1 0.3 2 0.4 3 0.6 4 0.8 5 06 1 0.2 0.1 0.2 0.2 0.2 0生成符合该分布随机数的步骤是：

1．生成一个（0，1）间均匀分布的随机数x0。

2．依据x0介于累积概率哪个区间来决定掷出骰子的点数x。如0

代码是

x0=rand; if x0<0.1 x=1;

elseif x0<0.3 x=2;

elseif x0<0.4 x=3

elseif x0<0.6 x=4

elseif x0<0.8 x=5 else x=6 end

这段语句能生成一个服从上表中离散分布的随机数x，如果生成多个x，可以用循环语句，也可以考虑将上述代码向量化。下图是我用上述代码生成10万个随机数所画出的分布直方图，可见这些随机数很符合上表中的分布律。 十七、生成多维联合分布随机数 一维随机变量是标量（也就是指单独的一个数字），而多维随机变量是一个向量。一个n维随机变量x是有n个分量的向量，(X0,X1,...,Xn)，用f(X0,X1,...,Xn)表示联合分布，用fk(Xk)表示第k维的边缘分布，用fk(Xk|X1=x1,X2=x2,....Xk-1=xk-1, Xk+1=xk+1,...,Xn=xn)表示当分量X1=x1,X2=x2,....Xk-1=xk-1, Xk+1=xk+1,...,Xn=xn时第k个分量xk的分布。这里大写X表示随机变量某个维度上的分量，小写x表示具体的数值。关于边缘分布、条件分布、联合分布

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库基于Matlab语言的Monte Carlo入门教程(6)在线全文阅读。