本文主要来自《全局光照技术》第五章。

蒙特卡洛方法是数值分析方法中的一个重要分支，它的核心概念是使用随机性来解决确定性的问题。大数定律告诉我们，满足某个概率分布的随机变量的数学期望所描述的积分函数可以使用这个随机变量随机抽样的样本均值来近似。

因为我们知道半径为1/2的面积是Pi/4，那么通过随机点的概率，就可以估算Pi/4的值，然后就可以估算出Pi的值。

随机变量X的每个值x都关联着一个概率，随机变量所有可能值组成的概率分布函数称为概率密度函数PDF，使用p表示。随机变量可以是离散的，对于离散的随机变量，每个随机数对应的概率满足：0 <= pi <= 1,所有离散变量的概率之和为1.

对于连续随机变量X，其概率密度函数p(x)是通过落于x附近区间[x,x+dx]内的随机数的概率p(x)dx来定义的，然而这种定义方式并不直观，所以离散随机变量的概率分布一般是通过累积分布函数CDF来定义，用大写的P来表示：

累积分布函数P(y)定义是所有随机数的值中小于或者等于y的随机变量的概率的积分，所以是一个递增函数，它具有以下属性：

随机变量的值x落于区间[a,b]的概率为：

对于离散随机变量X，假设其值为xi对应的抽样概率为pi，则该随机变量X的数学期望为：

期望代表的是对一个随机变量X进行抽样的平均结果。对应的，对于连续随机变量x，其期望为：

根据无意识的统计规律，设Y是随机变量X的函数，Y=g(X)，如果Σg(xi)pi绝对收敛，那么：

然后是方差的定义：

方差具有以下性质：

在统计学中，很多问题涉及对大量独立的随机变量的抽样xi的和进行处理，这些随机变量拥有相同的概率密度函数p，这样的随机变量称为独立同分布（IID）的随机变量，当这些随机变量抽样的和除以这些随机变量的数量N时，我们得到这些随机变量的期望的一个估计：

当N趋向于无穷的适合，我们可以确定随机变量统计的平均值趋近于期望的值：

一般最简单的求积分的方法便是将被积函数沿作用域上划分成多个区域，然后计算这些区域面积的和。如果这些区域的划分在作用域上是均匀的，那么对积分I的近似可以写为：

大数定律用于对期望表示的积分形式进行估计，即对积分函数进行估计，但是我们通常需要求得一个任意函数的积分，假设函数g(x)的定义域为x∈S，S可以是一个多维空间，此时x表示一个向量，我们希望计算以下积分函数：

而它的期望为：

对于上式进行一个变换，g(x)=f(x)p(x)，那么就变成：

我们定义I的估计In为：

In期望如下：

它的方差为：

所以随着N的增加，方差随着N线性降低。

蒙特卡洛方法可以归纳为下面的步骤：

考虑一个定义域空间Ω0以及一个概率密度函数p(x)，其中x∈Ω0，满足：

抽样过程是这样一个算法，它能够从p(x)对应的随机变量X中产生的一些系列随机数X1，X2...对于任意的Ω∈Ω0，满足：

由于我们只能通过随机函数求出一个均匀分布的值，那么怎么通过这些均匀分布的值产生符合概率密度函数p(x)的X呢？有以下几种办法。

逆变换算法

逆变换算法定义：设X是连续随机变量，其累积分布函数为Px,如果随机变量Y是一个[0,1]上的均匀分布，则随机变量Px-1(Y)具有和X一样的概率分布。

看这个图，dy表示随机变量是[0,1]上的均匀分布，那么p（x） = dP(x) / dx = dy / dx，Px的所以按照公式dP(P-1(Y)) / dx = dy/dx = p(x).

我们来看一个例子，一个离散的随机变量有4个可能值，概率为p1,p2,p3,p4.这些概率满足求和为1.

在一个标准的[0,1]上的均匀分布的随机变量分布在纵坐标上，通过在水平方向延伸可以和具有概率Pi的第i个输入事件相交，因为它是服从均匀分布，所以p3比p4更容易被相交，所以[0,1]上的均匀分布就被变换为服从概率密度函数p(x)的离散随机变量X。

所以步骤为：

取舍算法

考虑一个任意函数f(x)被限定在[a,b]区间，在区间外的所有值为0.

它的过程也很简单：

每次随机抽样产生的随机数矢量（X，Y）实际上是均匀的分布在一个矩形中，因此被接受的随机数矢量会均匀的分布与函数f(x)以下的区域。不过因为这样效率会比较低，所以我们通常使用要给建议分布来逼近f(x),

看一下新步骤：

原理也很好理解，本来是一个矩形，现在变成了一个三角形，自然被接受的概率就变大了。

马尔科夫链蒙特卡洛方法

考虑随机向量X~f(x1,x2.....xs)，每个分量的定义域为bi-ai，则直接抽样方法中取舍抽样的效率为：

其中L=maxf(x1,x2....xs)，如果L很大，或者s很多，那么拒绝率就很高。

为了解决这个问题，就有了新的算法。我们首先了解下马尔科夫链。设一个系统状态序列为x0,x1,x2....,每个状态可以以一定概率转移到另一个状态。假设我们从其中任意一个起始状态开始行走，每个时刻由一个状态转移到另一个状态，也可以转移到自身，对于任何时刻n，系统状态的转移条件概率为：

此转移过程形成的状态序列为X0,X1....Xn+1称为马尔可夫链。马尔可夫链的定义说明，将来时刻的n+1的状态，只依赖于当前时刻n的状态，与以前任何时刻的状态都无关。因此：

P(x0)为初始概率的状态，为1，P(Xn+1=x|Xn = y)称为转移率。如果转移概率和n无关，则成马尔科夫链是均匀的。通常我们研究的马尔科夫链都是均匀的。

马尔科夫链有几个属性.

1.不可约性。

如果一个系统从状态i开始，存在非零概率经过n>=1步之后可以达到状态j，则成i可以达到j。

表述为：

如果i->j， j->i同时存在，则称i和j是联通的，一个联通集合是指存在一个最大的状态集合C。里面每对状态都是互相连通的。如果一个状态空间就是一个单一的联通集合，则称该状态空间内从马尔可夫链是不可约的。

2.非周期性

假设从要给状态i出发，经过一定步数后回到自身，如果nii是k的倍数，则称状态i具有周期k，表述为：

这里gcd表示最大公约数，如果k=1，那么称为该状态是非周期态的，它意味着状态回到自身的步数可以是任意的，如果至少有一个状态是非周期态，则称为非周期的马尔科夫链。否则称为周期态。

3.回返性

假设从一个状态i触发，经过一定步数n回到i的概率记为：

概率Pii称为首达概率，如果平均返回事件Ti < ∞，那么称为状态i为正回返状态，否则，称为零回返状态。有限的马尔科夫链的回返状态必为正回返状态。

4.遍历性

如果状态i是正回返，不可约和非周期状态，则称状态是各态遍历状态。如果一个状态具有周期为1，并且能够在有限的平均步数内回返，则称为遍历态，如果有一个不可约的马尔科夫链中所有的状态都是遍历的，那么马尔科夫链是遍历的。

更一般的，如果要给马尔科夫链是遍历的，则始终存在一个步数n，使得从任意其他非i状态经过n步之后可以达到i状态，也就是说，马尔科夫链达到i状态的概率与初始状态无关。

5.平稳分布

设马尔可夫链在一个状态空间中从某个初始状态出发经过一定时间后，各个状态j的离散分布概率为Πj，如果对任意状态j存在概率密度函数Πj>0，

并且

即从所有其他状态转移到j状态的概率之和为状态j的概率分布：

则称马尔科夫链达到总体平衡。这里P是一个转移矩阵。一个正回返和非周期的马尔科夫链，当且仅当它是一个遍历链时，经过足够多的状态转移步数，这个马尔科夫链具有一个不变的概率分布Πj > 0, 这时无端初始分布如何，最终都会趋向于平稳概率分布Πj，于是有

我们用前面的天气为例，首先构造矩阵：

假设初始状态晴天为1，雨天为0.

经过一次随机行走后的状态概率分布：

经过n次运算，最终会这样：

达到稳定状态。

6细致平衡

在一个状态空间中，如果存在概率分布Π，使得：

该条件称为细致平衡，又称为局部平衡。细致平衡可以得到平稳分布。满足细致平衡的马尔科夫链都是可逆的，可逆性保证了平稳分布。

在这个过程中，最重要的时怎么根据已知的目标概率求出转移概率分布，梅特罗波利斯算法就是找出一种简单求解转移概率的方法。用马尔科夫链解决积分问题又两种不同思路：

1.已知目标概率分布pi，使用蒙特卡洛方法对这个目标概率分布进行抽样，然后计算积分。这种方法的关键时根据已知目标概率求出状态转移概率pij。

2.已知转移概率分布，然后给定任意一个初始概率分布，则马尔可夫链会根据转移概率逐渐趋向于稳定分布，这个一般用在辐射度方法中。

设f(x)是要给正比于目标概率分布p(x)的近似分布，梅特罗波利斯算法基本步骤如下：

那么怎么正面这个做法最终会收敛到唯一的稳定分布Π（x)呢？

首先我们先从满足细致平衡的条件开始：

梅特罗波利斯算法可以将转移概率计算分成两步，建议概率g(x'|x)和接受率α(x'|x)，使得：

根据前面的梅特罗波利斯算法可以知道它的选择为：

所以它是满足细致平衡条件的。

当然它也有自己的问题，我们为了加快收敛，就想要选择步幅大一点的建议密度函数，但这样会导致高频部分p(x')远大于p(x)，会让大量的随机数停留在尖锐区域。如果步幅过小，那么收敛就慢。

重要性采样

为了减少方差，我们可以使用重要性采样，它通过选择对一个与目标概率分布具有相似形状的分布函数进行抽样来减少方差。

一个完美估计的方差应该是0，因此我们希望能这样：

所以p(x) = |f(x)| / I. 然而我们不可能先算出I的值。我们可以选择和被积函数f(x)具有相似形状的概率密度函数来减少方差。直观上就是让f(x)/p(x)接近一个常数。

然后我们看下光照渲染中的直接光源计算：

我们可以选择f或者l来进行重要性采样。但结果会很差。

符合重要性抽样就是这样的一类抽样方法，他非常简单有效。他的步骤如下：

考虑如下的权重系数函数：

这里的ci是每个抽样分布pi对应的抽样数量的比例，ci = ni / N. ci是一个固定的数值，他在抽样之前确定。这个权重系数有一个属性，那么就是抽样值与i完全无关，因为每个x的抽样值对于所有的抽样技术都是无关的。从公式看其实就是pi(x)被约掉了。

不过重要性采样会导致有些区域被忽略。我们可以使用分层采样来处理这个问题。分层采样将被积函数的定义域Ω划分成彼此不相交的定义域Ω1...Ωn，对每个阶层内，固定数量的ni各随机数分别从抽样分布pi抽样而得。

为了加快蒙特卡洛方法的收敛速度，我们需要使用拟蒙特卡洛方法，他和传统的蒙特卡洛方法的不同在于随机数的产生方式。拟蒙特卡洛方法使用低偏差序列，这些序列是通过数论方法产生高度均匀的拟随机数序列。

首先我们要定义序列偏差度。

产生低偏差序列的方法和伪随机数序列的方法完全不同。具体算法就不展开了。

这次就到这啦。