《Computational Statistics with Matlab》硬译

第1章 从随机变量采样

　　研究者提出的概率模型对于分析方法来说通常比较复杂，研究者处理复杂概率模型时越来越依赖计算、数值方法，通过使用计算方法，研究者就不用对一些分析技术做一些不现实的假设（如正态性和独立性）。

　　大多数近似技术的关键是能够从分布中采样。需要采样来预测一个特别的模型在一些情景下是什么样的，找到在实验数据上应用模型的隐变量（参数）的合适值。大部分计算采样方法把从复杂分布采样的问题转化为简单采样分布的子问题。本章我们将介绍两种采样方法：逆变换方法和拒绝采样。这些方法适用于大多数单变量单值输出，下一章我们将讨论马尔科夫链蒙特卡罗方法，能够有效处理多变量分布。

1.1 标准分布

　　一些常用的分布就成了matlab支持的分布标准集的一部分。matlab统计工具箱支持大量的概率分布，使用matlab可以很容易计算概率密度、分布的累积密度以及从这些分布中采样，Table 1.1列举了matlab支持的一些标准分布，Matlab文档列举了更多可以用matlab模拟的分布。使用在线文档，很容易找到对其它常见分布的支持。

 

　　为了演示如何使用这些函数，Listing1.1展示了Matlab代码可视化Normal(µ,σ)分布，其中 µ = 100，σ = 15。为了帮助理解，想象这个分布代表人口的智商分布，代码演示了如何显示概率密度和累积密度，也演示了如何从分布中提取随机值，以及如何使用hist函数可视化这些随机采样的分布，代码产生的输出显示在Figure 1.1。

Listing1.1：可视化正态分布的Matlab代码

mu = 100;
sigma = 15;
xmin = 70;
xmax = 130; 
n = 100;
k = 10000;

x = linspace( xmin , xmax , n );
p = normpdf( x , mu , sigma );
c = normcdf( x , mu , sigma );

figure( 1 ); clf;

subplot( 1,3,1 );
plot(x,p,'k-')
xlabel( 'x' ); ylabel( 'pdf' );
title( 'Probability Density Function' );

subplot( 1,3,2 );
plot(x,c,'k-')
xlabel( 'x' ); ylabel( 'cdf' );
title( 'Cumulative Density Function' );

y = normrnd( mu , sigma , k , 1 );

subplot( 1,3,3 );
hist( y , 20 );
xlabel( 'x' ); ylabel( 'frequency' );
title( 'Histogram of random values' );

Figure 1.1 正态分布µ=100，σ=15

类似的，图1.2可视化离散Binomial(N,θ)分布，其中 N=10， θ=0.7。

代码，注意Binomial是离散分布，因此绘图用的是柱状图。

n = 10;
p = 0.7;
xmin = 0;
xmax = 10; 

k = linspace( xmin , xmax , 11 );
pdf = binopdf( k , n , p );
cdf = binocdf( k , n , p );
figure( 1 ); clf;
subplot( 1,3,1 );
bar(pdf,1,'c');
set(gca,'XTickLabel',{'0','1','2','3','4','5','6','7','8','9','10'},'FontSize',12)
axis([0 12 0.0 0.35])
xlabel( 'k' ); ylabel( 'pdf' );
title( 'Probability Density Function' );
subplot( 1,3,2 );
bar(cdf,1,'c');
set(gca,'XTickLabel',{'0','1','2','3','4','5','6','7','8','9','10'},'FontSize',12)
axis([0 12 0.0 1.0])
xlabel('k'); ylabel('cdf');
title( 'Cumulative Density Function' );
y = binornd( n,p,10000,1 );
subplot( 1,3,3 );
hist(y);
xlabel( 'k' ); ylabel( 'frequency' );
title( 'Histogram of random values' );

图1.2

练习：

1、采用listing1.1的matlab程序来展示Beta(a,b)分布，a=2，b=3。类似的显示指数分布Exponential(λ)，其中λ = 2。

a = 3;
b = 2;
xmin = 0;
xmax = 1; 
n = 100;
k = 10000;
x = linspace( xmin , xmax , n );
p = betapdf( x , a , b );
c = betacdf( x , a , b );
figure( 1 ); clf;
subplot( 1,3,1 );
plot(x,p,'k-')
xlabel( 'x' ); ylabel( 'pdf' );
title( 'Probability Density Function' );
subplot( 1,3,2 );
plot(x,c,'k-')
xlabel( 'x' ); ylabel( 'cdf' );
title( 'Cumulative Density Function' );
y =betarnd( a , b , k , 1 );
subplot( 1,3,3 );
hist( y , 20 );
xlabel( 'x' ); ylabel( 'frequency' );
title( 'Histogram of random values' );

mu = 2;
xmin = 0;
xmax = 5; 
n = 100;
k = 10000;
x = linspace( xmin , xmax , n );
p = exppdf( x , mu);
c = expcdf( x , mu);
figure( 1 ); clf;
subplot( 1,3,1 );
plot(x,p,'k-')
xlabel( 'x' ); ylabel( 'pdf' );
title( 'Probability Density Function' );
subplot( 1,3,2 );
plot(x,c,'k-')
xlabel( 'x' ); ylabel( 'cdf' );
title( 'Cumulative Density Function' );
y = exprnd( mu , k , 1 );
subplot( 1,3,3 );
hist( y , 20 );
xlabel( 'x' ); ylabel( 'frequency' );
title( 'Histogram of random values' );

2、用以上matlab程序来展示Binomial(N,θ)分布，其中N = 10，θ = 0.7。如图1.2。

3、写程序从Bernoulli(θ)中采样10个值，θ=0.3。Bernoulli分布是最简单的离散分布，只有两个输出值，0和1，以θ的概率输出为1，1-θ的概率输出为0。这个分布能模拟一些实验结果，如猜硬币的正反面，判断题的正确与否等等。在Matlab里，你可以使用Binomial分布来模拟Bernoulli分布，只需令N=1即可。尽管如此，为了练习的目的，请写出代码来模拟Bernoulli分布值，不要使用内置的Binomial分布。 

for n=1:10
randValue=rand(1)
if randValue<0.3
x(n)=1
else
x(n)=0
end
end

 4、保证每一次仿真都给出同样结果是非常有用的，在 Matlab 里，每次重执行代码，从分布中得出随机值是不同的。有一个简单的方法保证产生同样的序列，使用随机数字生成器。写从随机分布中采样10个随机值的Matlab代码时，在两次采样间使用seeding函数可使两个随机值集合相同。Matlab代码如下：

seed=1; rand(’state’,seed); randn(’state’,seed);

使用Matlab中的RandStream方法可以有更高级的方法来生成随机数。对于本文的大部分应用，用不到该方法。

5、假设我们从以前的认可研究中知道，智商系数是均值为100，标准差为15的正态分布。从该人口中随机抽出一人，计算他/她的智商大于110但小于130的概率。你可以使用一行Matlab代码。

 

normcdf(130,100,15)-normcdf(110,100,15)

 

**6、Matlab当前不支持Dirichlet分布，你能使用在线资源找到一个Matlab函数，来从Dirichlet分布中采样吗？

 1.1 从非标准分布中采样

      如果我们要从一个Matlab没有的非标准分布中采样，在建模时，这种情况经常遇到，譬如研究者能提出一个新的噪声过程或现有分布的组合。用计算方法来解决复杂采样问题，一般要依赖现有的简单分布采样方法。从这些简单分布中的采样值能被转换或同目标分布比较。其实本章讨论的一些技术已经在Matlab中用于复杂采样，如正态分布和指数分布。

1.2.1 离散变量的反变换采样

     反变换采样（也被称为反变换方法）是一种从概率分布中生成离散数字的方法，当然要给定累积分布函数的反。思想是从均匀分布随机数字（0和1之间）中采样，然后使用反累积分布函数转换这些值。这种过程的间接性是建立在基于转换均匀偏离的潜在采样的事实上。这种过程用于采样许多不同的分布。事实上，这就是Matlab执行许多分级数字生成器的方法。

    很容易将这种方法用在一个我们知道单个输出概率的离散分布。在这里，反转换方法只需要一个简单的查找表。

    给一个非标准离散分布的例子，我们使用一些实验数据，这个实验考察人类能够生成多均匀随机数字。在这些实验中，产生了大量的随机数字（从0到9），调查者将这些随机数字的相对频率制成表。正如你所想，人民并没有一直生成随机的分布。表1.2.1显示了一些典型数据。一些较小和较大的数代表性不足，而另一些数字则过多。由于一些原因，数字0和9从未产生（可能是对指令的解读有误）。这些数据说明人类并不擅长产生随机分布的数字。

    如果我们要模仿这个过程，写一个采样算法，根据表1.2中的概率。所以程序应该以0.2的概率生成4，以0.175的概率生成5，等等。例如Listing1.2使用Matlab内置的函数randsample来执行这个过程。代码生成效果如图1.2.1所示。

 

theta=[0.000;0.100;0.090;0.095;0.200;0.175;0.190;0.050;0.100;0.000]
seed=1;rand('state',seed);
K=10000;
digitset=0:9;
Y=randsample(digitset,K,true,theta);
figure(1);clf;
counts=hist(Y,digitset);
bar(digitset,counts,'k');
xlim([-0.5,9.5])
xlabel('Digit')
ylabel('Frequency');
title('Distribution of simulated draws of human digit generator');

 

注意上图在原文中是Figure1.3，但对该图的标题的错误的，明显不是BINOMIAL(N,θ)分布。第9页代码上面的Figure 1.2.1也是错误的，应为Figure1.4。

不用内置函数如randsample或者mnrnd，考虑如何使用反转换方法执行潜在采样函数是很有帮助的。我们首先计算累积概率分布。换句话说，我们需要知道一个结果等于或小于一个特定值得概率。如果F(X)表示累积分布函数，我们需要计算F(X=x)=p(X<=x)。对于离散分布，使用简单的求和就能完成。我们例子的累积概率显示在表1.2.1的右边一列。反转换算法的思想是从均匀离散偏离中采样，并同累积概率表中的每一个随机数比较，第一个输出的随机偏离小于采样输出相关累积概率。下图显示了一个例子，均匀随机偏离U=0.8，采样输出X=6。均匀重复采样偏离同累积分布比较，是离散变量反转换方法的基础。注意我们在使用逆函数，因为我们在做一个反查询表。

1、写一个Matlab程序对离散变量执行逆转换方法，用它来以表1.2.1的概率对随机数字采样。为了显示算法的运行，多采样些数字创建一个直方图。你的程序不能有0和9，因为它们在表中是0概率。代码如下，经测试，该代码生成的图形与上图一模一样。

cumpro=[0.000;0.100;0.190;0.285;0.485;0.660;0.850;0.900;1.000;1.000]
seed=1;rand('state',seed);
K=10000;
digitset=0:9;
for n=1:K
rd=rand();
for i=1:9
if rd>cumpro(i) & rd<=cumpro(i+1)
x(n)=digitset(i+1);
break;
end;
end;
end;

figure(1);clf;
counts=hist(x,digitset);
bar(digitset,counts,'k');
xlim([-0.5,9.5])
xlabel('Digit')
ylabel('Frequency');
title('Distribution of simulated draws of human digit generator');

2、上面练习的一个解决方法是通过使用多项式的随机数生成器 mnrnd，不需循环。如何使用这个函数来实现符合表1.2.1的采样。Matlab中未发现有该函数。

3、解释下为什么以上算法在处理陡峭的概率分布时效率很低。[提示：想象一种情景，前N-1个输出概率为0，最后一个输出概率为1]。你能找出一种简单的改变来提高效率吗？

1.2.2 连续变量的逆转换采样

逆转换采样方法也能用于连续分布。一般思想是均匀随机采样，然后将累积分布的逆函数用于随机偏离。下面用F(X)表示目标变量X的累积密度函数，F-1（X)是逆函数，假定我们可以计算这个逆。我们想要采样X值，可以通过以下步骤：

 

让我们用一个简单的例子来说明这个方法。假设我们要从指数分布中采样随机数。 λ > 0，累积密度函数是F(x|λ) = 1 − exp(−x/λ)，使用简单的代数，能够找到这个函数的逆，F−1 (u|λ) = −log(1−u)λ，这导致以下采样过程，从指数分数中采样随机数：

练习：

1、实现指数分布的逆转换采样方法，从这个分布中采样，并显示这些数值的分布。同指数分布的累积概率密度分布比较。

K=10000;
for n=1:K
rd=rand();
lamd=2;
x(n)=-log(1-rd)*lamd;
end;
hist(x,20);

**2、Matlab实现了一些函数，例如，当调用exprnd时，Matlab执行了存储在它自己内部路径中的函数。请定位Matlab函数exprnd并查看它的内容。Matlab如何从指数分布中采样？它是使用逆反方法吗？注意这个Matlab函数的路径依赖于你的Matlab安装，路径如下：

C:\Program Files\MATLAB\R2009B\toolbox\stats\exprnd.m

function r = exprnd(mu,varargin)
if nargin < 1
    error('stats:exprnd:TooFewInputs','Requires at least one input argument.');
end
[err, sizeOut] = statsizechk(1,mu,varargin{:});
if err > 0
    error('stats:exprnd:InputSizeMismatch','Size information is inconsistent.');
end
mu(mu < 0) = NaN;
r = -mu .* log(rand(sizeOut));

 可以看出exprnd正是采用逆反方法进行采样的。

1.2.3 拒绝采样

 很多情况下不能使用逆反方法采样，因为很难计算累积分布和它的逆。在这种情况下，有另外的选择，譬如拒绝采样，以及我们将在下一章讨论的马尔科夫链蒙特卡洛方法。拒绝采样方法的好处是它不需要任何"burn-in" period。所有的采样都能立即用于从目标分布中采样。

展示拒绝采样（也被称为接受拒绝采样）通用思想的一个方法是图1.5。假定我们希望得到圆心为(0,0)，半径为1的均匀点。首先，看起来从直接从圆中均匀采样非常复杂，尽管如此，我们可以用拒绝采样，首先从饶圆的正方形中获取(x,y)值，然后拒绝x^2+y^2>1的值得采样。注意在这个过程中，我们使用了一个简单的建议分布，如均匀分布，作为从更复杂分布中采样的基础。

      拒绝采样可以从很难采样但能评估一些特殊采样概率的分布中生成观测点。换句话说，假定有一个分布p(θ)，很难直接从这个分布中采样，但是我们能估计概率密度或者针对特殊值θ的概率p(θ)。研究者需要做的第一个选择proposal distribution，proposal distribution是一个简单的分布 q(θ)，能直接采样。思想是估计采用的样本在采用分布和目标分布下的概率，拒绝样本相对于建议分布，不可能在目标分布下。

     图 1.6 展示了整个过程。首先找到常量c，使对于所有的采样 θ，cq(θ) ≥p(θ)。常量c乘以建议分布 q(θ)作为比较分布，全部在目标分布上方。找到常量c并不困难，嘉定我们使用一些计算就能找到。现在从均匀分布[0,cq(θ)]中找出一个数字u，换句话说，这是0到我们提议的θ的比较分布的高度之间的一些点。如果u > p(θ)，我们将拒绝采样，其它就接受。如果我们接受采样，采样值 θ就是从目标分布 p(θ)中获得的值。这是这个计算过程的总结：

       这个算法的一个有高效操作的关键是包含尽可能多的采样。这严重依赖于建议分布的选择。一个建议分布与目标分布差别太大，将导致大量的采样被拒绝，减慢过程。

练习：

1、假定我们要从Beta(a,b)分布中采样，其中a=2，b=1。给定概率密度p(x)=2x，0<x<1。不适用matlab自带的Beta随机数生成器。目标是实现Matlab中的拒绝采样。作为练习，使用一个简单的均匀建议分布（即使不是一个好的建议分布）。常量c应为2，用直方图可视化采样值，检查分布是否契合Matlab自带的betarnd采样函数。接受的概率是多少？我们怎么改进拒绝采样器？

K=10000;
i=0;
x=[];
for n=1:K
randValue=rand(1);
rd=unifrnd(0,2*1);
if rd<2*randValue
x=[x,randValue];
end;
end;
hist(x,20);
length(x)

**2、图1.5的过程组成了Box-Muller方法的基础，生成高斯分布随机变量。首先生成均匀坐标(x,y)，从单位圆中，使用拒绝采样过程拒绝(x,y)对x^2+y^2>1。然后对于每一个(x,y)对，我们估计数量，z1和z2是每一个均为为0，方差为1的高斯分布。写一个Matlab程序，实现Box-Muller方法，验证采样值是高斯分布。

 

K=10000;
x=[];
for n=1:K
rv1=2* rand(1,1)-1;
rv2=2* rand(1,1)-1;
if (rv1^2+rv2^2)<=1
xxyy=rv1^2+rv2^2;
z1=rv1*(-2*log(xxyy)/xxyy)^(1/2);
z2=rv1*(-2*log(xxyy)/xxyy)^(1/2);
x=[x,z1];
x=[x,z2];
end;
end;
hist(x,30);

 第2章 马尔科夫链蒙特卡洛

数据概率模型的应用常需要用到复杂高维分布的积分，马尔科夫链蒙特卡洛是一种通用的计算方法，可以替代分析积分，通过累积从迭代算法生成的采样。用分析方法不能解决的问题，经常使用MCMC可以解决，即使是高维问题。MCMC的发展是统计中计算方法的最大进步。当MCMC成为一个非常活跃的研究领域，已经有了一些广泛使用的标准技术。本章我们将讨论两种MCMC， Metropolis-Hastings 和 Gibbs sampling ，在我们使用这两种技术之前，先了解MCMC的两个主要概念：蒙特卡洛积分和马尔科夫链。

2.1 蒙特卡洛积分

概率推断的许多问题需要在大亮结果空间的复杂积分或求和计算。例如，一个频繁问题是计算函数g(x)的期望，x是一元随机变量。如果x是连续纸，期望定义如下：

如果是离散变量，积分可用求和代替：

这些期望在许多情况下是上升的，当我们想计算一个分布的一些统计，如均值和方差。例如，g(x)=x，我们要计算分布的均值，使用分析技术积分或求和对于特定的分布就变得极具挑战性。例如，概率密度p(x)可能的函数形式不能使用分析积分。对于离散分布，结果空间太大，以至于对所有可能输出求和是不现实的。

蒙特卡洛积分是使用采样来估计一个复杂分布的期望。特别的，我们包含采样集合x(t),t=1,..,N，从分布p(x)中独立采样，在这种情况下，我们能使用有限和才估计2.1和2.2的期望：

在这个过程中，我们能使用合适的采样集来代替分析积分。一般情况下，估计的精度能够通过增加n来增加。估计的精度依赖于采样的独立性。当采样是相关的，采样大小增加的就有效。上一章讨论的拒绝采样不是个问题，但是一个与MCMC方法潜在的问题。

练习：

1、写Matlab代码来估计Beta(α,β) 的均值，使用蒙特卡洛积分，其中 α = 3， β = 4，你可以使用Matlab函数betarnd来从Beta分布中采样，你可以将你的答案跟分析答案 α/(α + β)比较。【提示：这只需一行Matlab代码】。

x=betarnd(4,3,10000,1);
mean(x(:))
3/7

2、类似的，估计Gamma(a,b)分布的方差，其中a=1.5，b=4，使用蒙特卡洛积分。Matlab命令gamrnd允许你从分布中采样。你的估计应该接近理论答案ab^2。

x=gamrnd(1.5,4,10000,1);
var(x(:))
1.5*16

2.2 马尔科夫链

 一个马尔科夫链是一个随机过程，从一个状态使用一个简单的序列过程转换到另一个状态。我们启动一个马尔科夫链，从一些状态x(1)，通过转换函数p(x(t)|x(t-1))，转到另一个状态，x(2)条件决定于上一个状态。我们通过迭代创建一个状态序列：

每一个这样的状态序列称为马尔科夫链。从马尔科夫链生成T个状态序列的过程如下：
1. Set t = 1

2. Generate a initial value u, and set x (t) = u
3. Repeat
t = t + 1
Sample a new value u from the transition function p(x (t) |x (t−1) )
Set x (t) = u
4. Until t = T
很重要的是，在这个迭代过程中，下一个状态只取决于上一个状态。所以，每一个马尔科夫链遍历状态空间，对下一个状态的转换只依赖于上一个状态。正是这种本地依赖造就了马尔科夫或无内存过程。我们将看到，这是马尔科夫链对于MCMC的一个重要特性。
当初始化每一个马尔科夫链时，状态空间将在起始状态周围。所以，如果我们启动许多链，每一个有不同的起始条件，这些链讲会接近于初始状态。这个时期称为收敛。马尔科夫链的一个重要特性是链的初始状态一旦经过充分长的转换序列后不再影响链的状态（假设马尔科夫链满足特定的给条件），在这点上，链达到了它的稳定状态，状态通过从稳定分布中采样反映出来。马尔科夫这种不管从什么地方开始都能收敛到稳定分布的特性是非常重要的。当应用MCMC时，它允许我们使用序列过程从分布中采样，起始状态也不会影响到估计过程。

例。图2.1是一个马尔科夫链的例子，涉及到一个连续变量。对于转换函数，从Beta分布采样。这个函数和它的常量使任意选择的，但是能够体现马尔科夫链的基本概念。这个过程从四个不同的初始值开始，每一个链持续迭代到100代。图显示了两个不同时间范围的序列。不同颜色的线代表四种不同的链。注意前10代显示了对初始状态的依赖，这是收敛过程，后面就稳定多了，如果我们不停止链将继续稳定。怎么能确定链收敛达到稳定状态了呢？尝尝不容易受，特别是高维状态空间。我们稍后将区分讨论收敛。

练习
1、开发matlab代码去执行例子中描述的马尔科夫链。创建一个类似于图2.1的面板。从0-1均匀分布中生成四个不同的初始值，启动马尔科夫链。[tip：……]

ix=unifrnd(0,1,1,5);
x=zeros(100,4);
for m=1:4
v=ix(m);
for n=1:100
x(n,m)=v;
v=betarnd(200*(0.9*v+0.05),200*(1-0.9*v-0.05));
end;
end;
plot(x)