MaxEnt模型

这里主要从物种分布建模（SDM）的视角出发，谈谈MaxEnt模型的原理和算法

引言

物种分布建模旨在利用已知的数据推断物种的在空间中的分布情况，这与最大熵原理存在相当大的共性。我们知道最大熵原理的核心就是：在所有满足已知约束条件的概率模型中，选择熵最大的那个模型。换句话说，当我们在对一个系统进行概率建模时，如果只掌握了部分知识（即约束条件），那么我们应该避免引入任何不必要的先验假设。最公平、最客观、最无偏的模型，就是那个在满足已知事实的前提下，内部不确定性最大的模型。该原理可以很自然的引入物种分布建模中，我们可以在已知物种出现情况和环境约束下去估计物种的分布。

数据

MaxEnt模型需要物种出现点数据以及背景数据。物种出现点数据即物种出现的位置数据，背景数据则是在研究区域中的采样，背景数据反映了整个研究区域的环境情况，是模型理解研究区域的基础，因此其非常重要。关于背景数据获取方法的内容不是本篇的重点，这里不多介绍，这里使用的背景数据采样方法为在研究区域的随机采样。

除了出现点数据和背景数据这些位置数据，该模型还需要环境变量数据，也就是每个位置点的各种环境变量值，包括但不限于：气温、湿度、土壤等等。需要指出的是，环境变量数据包括直接环境变量和间接环境变量，间接环境变量通过影响直接环境变量从而间接影响物种分布，因此在建模过程中最好使用直接环境变量。此外，要明确环境变量中的“环境”，比如生物相互作用不是其中的考虑因素；对于山脉、河流阻隔这些因素也不作为环境变量，它们被称为“Dispersal barriers”。

基本原理

MaxEnt模型即最大熵模型，基于最大熵原理，举一个例子：掷骰子

• 已知条件：你被告知一个骰子的平均点数是 4.5（即 $E[X] = 4.5$）。
• 问题：在仅知道这一条信息的情况下，如何估计骰子各个面朝上的概率 $p(1), p(2), ..., p(6)$？
• 错误做法：随意猜测一种概率分布，比如认为点数大的概率一定高。
• 最大熵做法：在所有满足 $\sum_{i=1}^{6} i \cdot p(i) = 4.5$ 和 $\sum_{i=1}^{6} p(i) = 1$ 的概率分布中，选择熵$H(p) = -\sum p(i)\log p(i)$最大的那一个，实际上就是一个带约束的优化问题。这个求得的分布，就是对我们未知情况最不冒险、最没有偏见的估计。

现在我们想要在一些已知条件下求解物种的分布情况，在引言中就提到了SDM建模与最大熵原理存在一些共性，我们能非常自然得到：

$$\begin{aligned} & \underset{p(x)}{\text{最大化}} & & H(p) = -\sum_{x \in L} p(x) \log p(x) \\ & \text{满足} & & \sum_{x \in L} p(x) = 1 \\ & & & \sum_{x \in L} p(x) h_j(z(x)) = \bar{h}_j \quad \text{对于 } j = 1, 2, ..., J \end{aligned}$$

其中：

• $H(p)$是分布的熵。
• $p(x)$表示已知物种出现时，这个出现发生在位置$x$的概率
• $\sum_{x} p(x) = 1$是概率的归一化约束。
• $h_j$是特征函数，$\bar{h}_j$是特征函数在物种出现点处的期望，$J$是特征函数的总量
• $z(x)$是位置$x$处的环境变量组合，例如（温度，降水，...）

这里点明几个可能会令人困惑的地方：

• 这里的自变量为位置$x$，在理论层面上$x$应该是研究区域内的任意位置，但是在实际操作上$x$是背景数据中的点
• $p(x)$的表示可能会有点拗口，让我们明确一下：$p(x)$的定义应该是$P(x|y=1)$，即在物种出现这个条件下，出现位置为$x$的概率，而我们真正想要求的是$P(y=1|x)$，即已知位置为$x$时，物种出现的概率。这二者通过贝叶斯定理联系$P(y=1|x) = P(x|y=1) * P(y=1) / P(x)$，其中的$P(y=1)$为物种普遍性，可以看作一个常量，$P(x)$为位置$x$被采样的概率，而每个位置被采样的概率是均等的（隐含假设，这是在认为背景无偏情况下的假设，上文提到了有不同的背景采样方法，在不同的方法下该假设未必成立）
• 这里的$\bar{h}_j$是依靠物种出现点的环境变量计算出来的
• 也就是说，背景数据提供了实际上的$x$的定义域，物种出现点数据提供了条件约束的特征值

求解结果： 使用拉格朗日乘数法求解这个优化问题，会得到一个非常优美且统一的解——指数族分布（也称为吉布斯分布或玻尔兹曼分布）：

$$p(x) = \frac{1}{Z(\beta_1, ..., \beta_J)} \exp\left( \sum_{j=1}^{J} \beta_j h_j(z(x)) \right)$$

其中：

• $\beta_k$是与第$k$个约束对应的拉格朗日乘子，需要通过数据来估计。
• $Z$是归一化常数（在统计物理中称为配分函数），确保所有概率之和为1。

上面是从地理空间出发，对于MaxEnt模型的理解，下面将从环境空间的视角来展开。

$$\begin{aligned} & \underset{f_1(z)}{\text{最小化}} & & D_{KL}(f_1 \| f) = \int f_1(z) \log \frac{f_1(z)}{f(z)} dz \\ & \text{满足} & & \int f_1(z) dz = 1 \\ & & & \int h_j(z) f_1(z) dz = \bar{h}_j \quad \text{对于 } j = 1, 2, ..., J \end{aligned}$$

其中：

• $D_{KL}(f_1 \|\| f)$是$f_1$和$f$的相对熵。
• $f_1(z)$是已知物种出现时，环境变量组合$z$的概率密度，$f(z)$是环境变量组合$z$在整个研究区域中的概率密度
• $\int f_1(z) dz = 1$是概率的归一化约束。

两种理解的数学表示实际上是一致的，两种表述等价的直观理解是：

• 在地理空间中，最大熵分布相对于均匀分布
• 在环境空间中，这个分布对应的$f_1(z)$就是与$f(z)$相对熵最小的分布

数学上，可以通过以下步骤证明：

1. 地理空间的均匀分布$u(x) = 1/|L|$对应于环境空间的分布$f(z)$
2. 地理空间的任意分布$p(x)$对应于环境空间的某个分布$f_1(z)$
3. 可以证明$D_{KL}(p \| u) = D_{KL}(f_1 \| f) + \text{常数}$
4. 因此，最小化$D_{KL}(p \| u)$等价于最小化$D_{KL}(f_1 \| f)$

在环境空间表述下，我们想要求的是条件概率$p(y=1|z)$，即在给定环境$z$下物种出现的概率，根据贝叶斯公式有：

$$p(y=1|z) = \frac{f_1(z) * p(y=1)}{f(z)}$$

其中：

• $p(y=1)$是物种普遍性，在仅出现数据下是不可辨识的（即无法仅从出现数据中得知）。
• 因此，我们无法直接得到$p(y=1|z)$，但可以估计$f₁(z)/f(z)$。

解惑

这里总结一些在学习过程中的疑惑。

对于一些变量的明确

首先，我们求解得到的$p(x)$不是物种在位置$x$出现的概率；而是已知物种出现的条件下，该出现位置为$x$的概率。这里可能会有疑惑，为什么不直接求物种在位置$x$出现的概率呢？这时由模型使用的数据所决定的。

模型使用了这些数据：

• 出现点数据：$\{x_1, x_2, ...,x_n\}$，这些都是已知$y=1$（即物种出现）时的$x$值
• 这些出现点数据正式从$P(x|y=1)$的分布中抽取的样本

我们没有：

• 在不同位置$x$处物种出现/不出现的观测数据

其次，$p(x)$，即$P(x|y=1)$无法直接计算得到真正的$P(y=1|x)$，在前面我们提到了：

$$P(y=1|x) = \frac {P(x|y=1) * P(y=1)} { P(x) }$$

对于$P(x)$，理论上它只依赖于研究区域$L$的规模，在认为背景数据无偏时，我们通过背景数据反映整个研究区域，则$P(x)$的值只依赖于背景数据的规模，当背景数据确定时，其值也就确定了。而$P(y=1)$为物种普遍性，这是一个常量，但是是一个未知常量，因此在该值未知的情况下，我们无法通过$P(x|y=1)$直接求解出真正的$P(y=1|x)$分布。也即我们无法求出物种在位置$x$出现的绝对概率，但是，我们可以求解出一个$P'(y=1|x)$分布，且该分布与$P(y=1|x)$具有相同的适宜性。也就是说，对于任意两个位置$x_1,x_2$， 我无法知晓这物种在这两个位置出现的具体概率$P(y=1|x_1)$和$P(y=1|x_2)$分别是多少，但是我能够确定$P'(y=1|x_1) > P'(y=1|x_2)$，那么就有$P(y=1|x_1) > P(y=1|x_2)$。这对于大多数应用是足够的，我们只需要知道哪里比哪里更好，而不需要知道绝对的出现概率值。这一点在下面从环境空间视角出发也是同样的。

关于$f(z)$和$f_1(z)$

$f(z)$表示背景分布，是对于整个研究区域环境的反映。如果不考虑计算量、边际效应等因素，我们对整个研究区域进行采样，也即背景数据包含了完整的研究区域，这样得到的$f(z)$是最无偏、最准确、最能真实反映研究区域环境分布的。但是在实际操作层面上，我们不可能对于整个研究区域做完整的采样，一般都是对研究区域做随机的采样，以采样到的有限样点来估计$f(z)$。事实上，如果是无偏的随机采样，一定数量的样点构成的背景数据从而估计出的$f(z)$已经能比较好的反映整个研究区域的环境分布了。

经验分布误区

$f(z)$描述了研究区域提供的环境，而$f_1(z)$则描述了物种所利用的环境。这里可能会存在一些疑问，我们直接通过背景数据能够统计计算得到$f(z)$,那么直接对于物种出现点数据做统计计算不就得到$f_1(z)$了吗？

事实上，这是问题的起点，而不是终点。$f(z)$通过统计得到的是经验分布，当背景采样无偏且充分时，根据大数定律我们知道该经验分布会接近于真实分布，也即$f(z)$能较好地反映研究区域的环境分布情况；类似的，如果我们通过统计得到$f_1(z)$的经验分布， 我们会发现该分布完全满足任意环境约束，但是它与$f(z)$的相对熵未必最小。

此外，经验分布看似合理，它其中还包含了若干问题：

• 我们不能直接以经验分布作为结果，它没有考虑基准。举一个比较极端的例子，比如我统计了某个学生每天午餐的情况，我发现他每天的午餐点的全是青菜，所以我得出了该学生偏爱青菜的结论，但是该结论是没有考虑基准的，我应该还统计每天中文食堂提供的饭菜，如果食堂每天中文提供的只有青菜，那么学生也就只能点青菜，那么之前得出的偏爱青菜的结论就是有问题的。
• 经验分布还可能存在过拟合和泛化能力差的问题，它可能无法预测为观测的数据

最大化熵和最小化相对熵

在地理空间表述中，我们的目标是在给定环境约束下，求出一个物种在空间中的熵最大的概率分布$p(x)$；而在环境空间表述中，我们的目标是在给定环境约束下，求出与背景分布$f(z)$最接近的$f_1(z)$，即最小化$f_1(z)$与$f(z)$的相对熵。

这两种表述是等价的。实际上它们的模式其实是一致的，在地理表述中我们求的是熵尽可能大（要满足约束）的概率分布，而我们知道熵最大的的分布就是均匀分布，实质上我们也是要使得我们求的分布尽可能接近均匀分布，即最小化$p(x)$与均匀分布$u(x)$的相对熵，但是一般不这么表述，而是直接称最大化熵。

而在地理空间中的均匀分布$u(x)$与环境空间的的背景分布$f(z)$是相对应的。需要指出的是，这里的$f(z)$理论上描述了整个研究区域的环境分布，但是在实际上它是通过背景数据计算得来的，但是在这里我们只需知道其在理论上的表征即可，它描述了环境$z$在整个环境空间中出现的概率为$f(z)$。不难看出，如果你认为物种在地理空间中是均匀分布的，那么当你将地理位置向环境空间中映射时，就会发现物种在某个确定的环境$z$中出现的概率，正比于该环境在整个研究区域中出现的概率（即$f(z)$）。

在啰嗦一些：

在地理空间中，我们要求$p(x)$，即已知物种出现的条件下，该出现位置在$x$的概率，我们要使得该分布尽可能接近均匀分布$u(x)$；

在环境空间中，我们要求$f_1(z)$，即已知物种出现的条件下，该出现环境为$z$的概率，我们要使得该分布尽可能接近研究区域的环境分布$f(z)$。

这样来看这二者的模式就非常一致了

背景数据

背景数据的选取是MaxEnt的基石，$f(z)$由背景数据得来，而$f(z)$又是模型的比较基准。改变了基准，就改变了“正常”的定义，从而彻底改变了模型对“偏好”的判断。需要指出的是，背景数据是研究区域的反映，在这一段文字中，这里的背景数据与研究区域可以看作等价的，即研究区域规模的选择对于模型结果的影响也是巨大的，例如在论文（A statistical explanation of MaxEnt for ecologists）的研究案例中有：

• 当背景是整个澳大利亚时，模型认为 等温性 是最重要的变量，因为它能很好地将物种分布的两南角与澳洲内陆的广袤干旱区区分开。
• 当背景是局域的西南植物区时，模型则认为最干季降水和最热季温度更重要，因为这些变量能解释物种在局域范围内的分布差异。

有偏背景

背景数据对于MaxEnt模型至关重要，因此这里额外分出一些章节来说明。

如果我们认为出现点是有偏的，那么我们应该使用有偏的背景数据，且这二者的偏差应该是尽可能一致。

其原理在于： MaxEnt（以及其他仅出现模型）的核心是比较$f_1(z)$（出现点的环境分布）和$f(z)$（背景点的环境分布）。如果出现点数据因为采样偏差（例如，只沿着公路采集）而扭曲，那么$f_1(z)$估计的其实是$f_1(z) * s(z)$，其中$s(z)$是采样偏差函数。

• 如果使用无偏背景：模型会错误地将采样偏差$s(z)$解释为物种的生态偏好。模型会认为：“物种更喜欢公路附近的环境”，而实际上可能只是“调查者更喜欢在公路附近调查”。
• 如果使用有偏背景：让背景数据来自同样的偏差采样过程（即背景分布也乘以同一个$s(z)$）。那么，模型比较的是$[f_1(z)s(z)] / [f(z)s(z)]$，偏差项$s(z)$在分子分母中被抵消了。模型最终学到的是真正的$f_1(z)/f(z)$，即物种的真实环境偏好。

结论：当出现点数据存在已知或可估计的采样偏差时，使用具有相似偏差的背景数据是一种有效的统计校正手段。

上面是从环境空间的视角出发考虑背景偏差的，从地理空间上也是类似的，在一般情况下，我们要使得$p(x)$尽可能接近均匀分布$u(x)$（即熵最大化），然而在认为出现点有偏的情况下，估计的$p(x)$其实是$p(x)*s(x)$，因此我们需要对背景做相同的处理，即$u(x)*s(x)$，也即有偏的背景采样，这里就不再是随机采样了。

特征函数与正则化

特征函数，即$h_j(z)$是约束条件的重要组成，而正则化是防止模型过拟合的方法，在这里正则化方法影响了特征函数的系数。

对于特征函数，MaxEnt通过预先生成大量候选特征，然后在优化过程中自动选择重要特征来解决这个问题。

第一步：特征生成（预处理）

在模型拟合之前，MaxEnt会为每个环境变量预先生成大量的候选特征：

• 对于连续变量：
  • 生成多个铰链特征，铰链点设在数据的各种分位数处
  • 生成多个阈值特征，阈值设在数据的各种分位数处
  • 生成线性特征和二次特征
• 示例：对于温度变量，可能生成：
  • 铰链特征：max(0, temp-10), max(0, temp-15), max(0, temp-20), ...
  • 阈值特征：I(temp>10), I(temp>15), I(temp>20), ...
  • 线性特征：temp
  • 二次特征：temp²

第二步：约束的重新表述

现在，约束条件不再是关于"带参数的特征"，而是关于这些预先定义好的具体特征：

我们需要找到$f_1(z)$，使得对于每一个预先生成的特征$h_j$：$E_{f₁}[h_j(z)] = \bar{h}_j$

第三步：通过正则化进行特征选择

这是最关键的一步！MaxEnt使用**L1正则化（套索）**来自动选择哪些特征应该被包含在最终模型中，这里的自动选择实际上就是通过$h_j$的系数$\beta_j$来控制。

综上，我们发现实际上特征函数的数量非常的多，但是最后只会有部分特征被保留到模型中。

对于正则化，这里可能涉及比较具体的运算层面了。我们知道正则化通过引入惩罚项来防止模型过拟合，在MaxEnt模型中，我们一般选择L1正则化，因此损失函数发生了变化，我们原本是要最小化$f_1$与$f$的相对熵，而现在变成了最小化：

$$L=D_{KL}(f1 || f) + Σ λ_j |β_j|$$

其中：

• $Σ λ_j |β_j|$是惩罚项
• $\beta _ j$是特征函数$h_j$的系数，事实上$\beta_j$是每个约束所对应的拉格朗日乘子
• $λ_j$是每个特征对应的正则化强度参数，一般来说正则化参数越大，对应的系数参数越小
• 这里有$λ_j = λ × \sqrt {s²[h_j] / m}$其中的$λ$为超参数，$s²[h_j]$为$h_j$在出现点上的方差，$m$为出现点数据的数量

上面展示了在地理空间中MaxEnt模型的解，通过转换，在环境空间中的解形如：$f_1(z)=f(z)exp(\sum_j\beta_j h_j(z))/Z$。不难发现我们需要的$f_1$是依赖于$\beta_j$的，也就是这些拉格朗日乘子，事实上该优化问题的解就是这些乘子，需要通数值方法来求解。

迭代的过程大致为：

a. 根据当前 β，计算$f_1(z) = f(z)exp(∑β_j h_j(z))/Z$ b. 计算每个特征的模型期望$E_{f_1}[h_j]$ c. 对于每个特征$j$：

• 如果$|E_{f_1}[h_j] - \barℎ_j| < λ_j$：设置$β_j = 0$
• 如果$|E_{f_1}[h_j] - \barℎ_j| = λ_j$：保持当前$\beta_j$的值
• 否则：增大$β_j$以减少$|E_{f_1}[h_j] - \barℎ_j|$

对于上述过程：

• 我们知道$f_1(z)=f(z)exp(\sum_j\beta_j h_j(z))/Z$，当我们增大$\beta_j$时，实际上是增大了特征$h_j$在$f_1$中的比重，这强化了模型对于对应特征的捕捉能力。因此$|E_{f_1}[h_j] - \barℎ_j|$会减小，当$|E_{f_1}[h_j] - \barℎ_j| = 0$时表示模型完全捕捉了给定数据的对应特征，然而为了防止过拟合我们引入了惩罚项，为平衡正则化惩罚我们一般设置最优停止点为$|E_{f_1}[h_j] - \barℎ_j| = \lambda_j$。
• 对于$|E_{f_1}[h_j] - \barℎ_j|$和$\lambda_j$的关联源于损失函数$L$关于$\beta_j$的梯度：$\frac{\partial L}{\partial \beta_j} \propto (E_{f_1}[h_j]-\bar h_j)+\lambda_j \cdot sign(\beta_j)$，该梯度反映了当$\beta_j$变化时损失函数$L$变化的速率，这个梯度包含两个部分，第一部分为相对熵的梯度（收益），后一部分为惩罚项的（成本）。该式子表明当我们调整$\beta_j$时，我们将会引入$|E_{f_1}[h_j] - \barℎ_j|$的相对熵减少（也就是收益，毕竟我们希望减小损失函数），而会得到$\lambda_j$增大的惩罚（也就是成本），因此我们需要权衡二者，确保收益大于成本时进行调整。（这里的增大不是增大的绝对的值，而是一个相对的比例，这里举一个简单的例子来辅助理解，对于一个简单的函数$g(x)=ax$，有$g'(x)=a$，这说明当$x$变化$\Delta x$，总会获得值为$a\Delta x$的增量，我们可以控制$\Delta x$来控制增量，但是$a$规定了增量的比例，这里的收益和成本本质也是一种比例）