蓝色的图形就是用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表示具体的数值。关于边缘分布、条件分布、联合分布一定要明白，这些都是基础数学知识，非本课程内容。如果手头没有书，通过google搜索或上维基百科临阵磨枪也是可以的。 各种生成多维分布随机数的方法一般步骤都是，逐个维度生成随机数分量，最后将这些分量依次组合起来——如先生成x0,再x1,...,最后xn,，最终写成(x0,x1,...,xn)。 在详细讲如何生成这些分量前，我们讲讲如何储存生成的随机数。 如果一次生成一个n维的随机数向量，可以用n变量来储存这个随机数的n个分量，也可以将这n个分量按照次序（次序不能乱）存于一个1*n的行向量中。如果一次生成随机数的数量很多，例如N个随机数，前面两种办法都可以用，即可用n个变量来储存这些随机数的每个分量，此时每个变量是N*1的列向量；也可以只用一个N*n矩阵储存随机数所有 分量，这个矩阵每一行是一个服从规定的联合分布的随机数，共有N行即表示共储存N个这样的随机数，矩阵的每一列表示这N个随机数的一个维度上的分量，共有n个维度。 （一）最简单的——各维度独立 各维度独立的联合分布随机数的生成最为方便，由于联合分布函数就是每个维度边缘分布函数的直接乘积，所以只要分别生成每个维度的随机数分量，然后组合成随机数向量即可得到服从该联合分布的随机数。 例子1，生成一个在0≤x≤2,0≤y≤2，正方形区域上的二维均匀分布。二维均匀分布在每个维度上都是均匀分布（即两个维度的边缘分布都是(0,2)上的均匀分布），且两个维度互相独立。 用第一种存储方法， x=unifrnd(0,2); y=unifrnd(0,2); 则每个维度上分别生成一个服从(0,2)均匀分布并分别储存在x,y这两个变量中。如果一次生成多个随机数，如N个，可用 N=400; x=unifrnd(0,2,N,1); y=unifrnd(0,2,N,1); 这里x，y都是N*1大小的列向量，分布存储着这N个随机数的第一维和第二维两个分量。我们看看这些随机数是否很好得符合二维均匀分布特性。 scatter(x,y); 接着，我们看看用上述第二种存储方法， X=[x,y]; 紧接着第一种存储方法中的语句，我们将生成的两个分量组合起来储存到一个变量中。当然这里还有一种取巧的办法，由于两个维度的边缘分布都相同且独立，我们只需用unifrnd函数一次性生成一个N*n大小的矩阵就可以了。 X=unifrnd(0,2,N,2); 例子2，我们要生成的随机数服从一个三维联合分布，其第一维边缘分布服从标准正太，第二维边缘分布是自由度为4的t分布，第三维边缘分布是自由度为(7,8)的F分布，各个维度边际的边缘分布之间相互独立。我们只要用如下代码： x1=rand x2=trnd(4) x3=frnd(7,8) x=[x1,x2,x3] x1,x2,x3分布储存三个维度的分量（第一种存储方法），将这些分量组合起来存入x中（第二种存储方法）。 如果要一次就能生成一堆这样的随机数。可以用如下的代码： N=1000; x=[rand(N,1),trnd(4,[N,1]),frnd(7,8,[N,1])]; 这段代码略过了中间过程，直接生成第二种存储方法所说的矩阵，这个矩阵大小为N*3， 我们可以大致观察该联合分布在每个区域内的概率密度的大小。 scatter3(x(:,1),x(:,2),x(:,3),'marker','.','sizedata',1); 注：x(:,1)表示将x矩阵的第一列（也即随机数第一维上的分量）提取出来。
蒙特卡罗算法与matlab(精品教程)
