BME（Bayesian maximum entropy）

最近在看BME，看了很多文章，有这样一种感觉，就是理论都看懂了，但是不知道咋使用，不知道在实际运用过程中是咋样一步步操作的，在此记录下BME的理论和实践的学习过程，方便自己梳理。

1.BME理论基础

1.1 一些变量说明

首先，大自然的时空发生过程可以看做是一个随机过程（The spatiotemporal random field ，STRF），然后在一些时空点位(p1,p2,…,pn)上，有着对应的随机变量（x1,x2,…xn)，这些随机变量产生了一些实例（X1,X2,…,Xn），其中pi=(si,ti)，si是n为空间位置，ti是时间，这就是大自然产生数据的过程。也就是BME的一个思想框架，BME就认为这些大自然产生的这些数据是由这些时空点位上的随机变量生成的，因此满足一定的概率密度函数，可以用后验概率分布来描述它们，以此来进行预测。

1.2 数据说明

这个框架中广义的说数据分成3种，第一种是知识，知识可以是物理定律，可以是人们的经验，可以是其他的知识，都可以；第二种是hard data，hard data就是可以精确描述的数据，比如人的身高数据，1.4m,1.6m之类的这种，其他的模型，如机器学习模型，就是使用的hard data；第三类是soft data，这类数据就是和hard data相对应的，不能精确描述的数据，经典的soft data包括区间数据、分布数据、函数数据，比如我们知道某个人的身高是在[1.5,1.8]m之间，这种就是区间数据。

1.3BME框架/流程

1.3.1 Prior stage

先验阶段实际上是借助硬数据和知识来建立先验概率密度函数（pdf），建立先验概率密度函数的过程实际上求解一个最优化问题。
目标函数
目标函数是信息熵，希望能使信息熵最大，也就是包含的信息尽量多，在BME中就是希望pdf能够尽可能多的包含数据的信息，理论上只有目标函数的话均匀分布是使得信息熵最大的。

约束
在上面的过程中似乎没有任何数据和知识，按照熵最大的目标得到了均匀函数。但是均匀函数可能和我们的知识相当不吻合，比如我们的知识告诉我们这个分布是高斯的，比如我们的知识告诉我们某个区间的概率是0.3，比如我们的知识告诉我们小于某个阈值的概率是0.7，这些知识就可以成为约束条件。
除了知识外，硬数据的信息也会成为约束，具体而言，硬数据的各阶矩一般都需要满足，比如一阶矩（均值），二阶矩（方差），以及高阶矩等等，通过硬数据我们可以算出样本一二三等阶矩，那么我们要求广义分布要和这些样本（硬数据）相一致，具体用几阶矩可以视情况而定。
于是目标方程和约束就构成了一个优化问题，使用拉格朗日乘子法可以求解这个优化问题，得到先验函数的通式如下所示：

其中的$u_0$和$u_{\alpha }$是未知的系数，在后验概率分布求解阶段求取。

1.3.2 Metaprior stage

这个阶段主要是对数据做处理，包括数据清洗、转换等等，准备hard data和soft data的形式，比如将soft data转化为之前提到的三种形式之一（区间软数据、概率分布软数据、函数软数据），为后验概率分布求取做好准备工作。

1.3.3 Posterior stage

后验概率分布求解阶段就是要根据软硬数据去求取后验概率分布，后验概率分布求取的通式如下所示，这个其实就是一个条件概率公式，已知data的条件下，知道各个位置处的分布。

对于区间软数据，后验概率分布计算公式如下：

对于概率软数据，后验概率分布计算公式如下，式中的F（.）是累计概率函数。

---

目前我的理解是，Posterior stage阶段是利用软硬数据求解出$f_G({\chi }_{map})$中的未知的系数，类似于机器学习模型中的训练阶段，不过这里有点特殊的是训练不仅仅用到了hard data，还用到了soft data。硬数据的训练方式我的理解是将自变量代入后得到了一个因变量的分布，然后因变量分布的均值或中位数或众数等于硬数据中的因变量，以此去训练参数。软数据的训练方式我的理解是代入自变量后（以区间软数据为例），得到一个因变量的分布，然后对这个分布进行软数据区间积分，出现在这个区间的概率等于置信区间，以此来训练获得参数。
但是目前我仅仅是看了理论，还没有具体实践过，因此对这个理解不敢确定，很有可能是错的，若大佬知道可以指出来。

2.基于贝叶斯最大熵的多因子空间属性预测新方法

看了一篇具体的文章，看完了还是不知道细节是咋计算的，也找不到资料，好烦恼。
文章使用的是基于克里金方法和环境相关法的BME模型，该模型可以比单纯的使用克里金法效果更好。
下面讲一下在这篇文章中是如何去使用BME方法的。

2.1先验阶段

首先可以证明根据贝叶斯最大熵原理,当 $K_G$仅考虑硬数据的各阶统计动差时,所得结果与普通克里金法的结果一样;因此,在本研究中先验 pdf 可用普通克里金方法计算, 结果符合正态分布, 即表示为:

根据这篇文章的说法，先验应该是可以完全求解出一个分布的，并不存在未知参数，这个例子里的先验函数并不存在任何未知数，这样看来的话我上面的理解就错误了。这个地方还是先mark一下，究竟如何还需要探究。

2.2中间阶段

中间阶段是准备软数据的阶段，这边的软数据使用的是环境相关法得到的，下面简单介绍一下环境相关法，这个方法也相对简单。
假设有一个自变量EA，因变量是CA，分别将自变量和因变量分为N和n组，计算出类别-环境因子定量关系为：
$$[(\frac{Count(1{)}_i}{Coun{t}_i}，C{A}_1),...,(\frac{Count(N{)}_i}{Coun{t}_i}，C{A}_n)]$$
其中 $Coun{t}_i$为所有建模样点上辅助变量 EA的取值属于 $E_i$ 的个数，而$Count(k{)}_i$ 为 $C{A}_K$ 包含的样点上辅助变量 EA 的取值属于 $E{A}_i$ 的个数。
上面这个计算出来的实际就是自变量取值是$E{A}_i$这个区间的时候，因变量CA的分布是咋样，在这篇论文里这个描述似乎不大准确。这样一来通过自变量我就可以知道因变量的分布是咋样的了。但是这是一个自变量的情况，多个自变量就是对单个自变量做加权，权重使用的是自变量和因变量的相关系数，其中注意需要做归一化，因为最后分布的概率密度之和是1.
通过上述的环境相关法就得到了一些未测点的需要预测点的概率分布，这个也就是软数据。

2.3 后验阶段

后验阶段依然是我最不能理解的一个阶段，这个阶段没有找到任何有细节描述的资料，论文也只是给出了一个general的条件概率公式：

但是我还是不知道到底细节是咋操作的，如何使用软硬数据去更新概率密度函数，得到后验概率分布。

3.软件实现

可以通过实践过程来更进一步的了解BME模型，软件实现使用的是BMEGUI3.0.1，软件地址点击该链接。里面包含软件安装所需要的library，该软件包，还有使用手册、安装手册、示例数据集，各示例数据集的tutorial。

3.1软件安装

下面记录下软件安装过程中所踩的坑，实在有点头秃：

1. BMEGUI软件在不断更新，可以看一下现在更新到第几版了，下一个最新版的，我给出的链接是3.0.1版的，可以点击这里自己去看看有没有更新的版本。
2. 按照安装手册说明的顺序先将包一个个的安装好，注意安装过程中退出杀毒软件，比如360，否则过程中可能有bug，比如我在安装MCR7.13的时候就因为杀毒软件安装失败了，报的错还查不到原因，折腾了好久。注意一定要每个包都安装好了。
3. MCR安装好后，运行BMEGUI报错说找不到对应版本的MCR，检查该版本是否安装成功，然后检查系统路径中有没有MCR，如果都没问题，但是依然报错，重启电脑试试。
4. 下载下来的BMEGUI安装包不要放在有中文路径的位置下，保持路径全英文，要不然会出bug。
5. 注意电脑上可能同时装了python2和python3，而且两个都在系统路径下，此时很可能默认是python3，这时需要做一下调整，我的方法是找到python2的命令所在的文件，然后将python.exe修改成python2.exe，这样一来cmd下面打python2就是调用的python2了，和python3就区分开来了。
6. 按照安装手册说的双击.py文件，可能直接打开了源代码，而并没有执行，这时右键，选择打开方式，在打开方式中找到python2.exe所在文件夹，然后点击python2.exe，此时.py文件的默认打开方式就变成python2了。
7. 双击00 create_desktop_GUI.py文件就可以在桌面创建图标了，双击图标就可以执行了。
8. 具体使用过程中，注意工作目录路径、名称和文件都不要出现中文，所有的都一样是英文，要不然直接无法进入下一步，连错误提示都不会有，文件格式严格按照要求，要不然也是直接无法进入下一步，连错误提示都不会有。

3.2软件输入数据

数据包括一个是分析数据源一个是river network data，其中分析数据源包括时间空间和观察数据，它既可以是软数据也可以是硬数据，其中硬数据比较容易，软数据可以包括区间数据、高斯分布、三角分布、截断高斯分布，当同时包括软硬数据时，需要加一个type字段来说明是什么类型的数据。river network data直接看下面的示意图就可以很好理解了。

需要说明的是为什么要使用river network data，因为度量距离的时候默认是用经纬度计算直线坐标，但这很多时候不符合实际情况，加上河网拓扑之后，就是使用拓扑来计算距离了，两点间网络最小距离作为两点间的距离，不过BMEGUI并不是使用图结构来刻画的，而是使用的有向树，而且必须指明outlet，因此之后如果想要勾画路网，可能还要想想到底适不适用。
具体说明可以参考用户手册了。

3.3 Dialog Box 2 (Data Distribution)

主要就是一些描述性统计的结果，直方图，均值、方差等统计指标了，这个没啥好说的，也非常容易理解。
可以说一下的是，由于我们处理的数据中还包括软数据，软数据是没法做描述性统计之类的，这个实际上是先对软数据进行硬化，然后使用“Hardened” values去进行描述性统计之类的，硬化方式就是取均值。

3.4 Dialog Box 3 (Exploratory Data Analysis)

只是画个图看看数据咋样的而已，没啥好说的，直接贴手册原话了

Dialog Box 3 (Exploratory Data Analysis) is used to conduct the exploratory data analysis. This dialog box has two tabs, labeled “Temporal Evolution” and “Spatial Distribution”, respectively. BMEGUI displays the time series plot of the measurement values at each monitoring location on the “Temporal Evolution” tab, and the spatial distribution plot of the measurement values at specific times on the “Spatial Distribution” tab.

然后用户可以用Aggregation Period去对时间做聚合，然后展示，也是简单的。

3.5 Dialog Box 4 (Mean Trend Analysis)

全局均值趋势分析使用的下面的公式：

其中$m_{ss}(s)$代表的是空间分量，$m_{ts}(t)$代表的时间分量。

1. 对每个观测点取均值得到the raw spatial mean $m_s$;
2. 使用指数空间滤波器得到spatial component $m_{ss}$ that is smoothed over space；
3. 对每个观测时间取均值得到the raw temporal mean $m_t$;
4. 使用指数时间绿波器得到temporal component $m_{ts}$ that is smoothed over time。

这里简单说一下exponential filter，exponential filter公式如下图所示：（Exponential smoothing is a rule of thumb technique for smoothing time series data using the exponential window function. Whereas in the simple moving average the past observations are weighted equally, exponential functions are used to assign exponentially decreasing weights over time. ）

具体可见wiki。
在BMEGUI中，需要输入的参数包括4个：“Spatial Search Radius”,“Spatial Smoothing Range”，“Temporal Search Radius”, “Temporal Smoothing Range”，我理解的这些参数的意思是：首先Temporal Search Radius的意思比较容易理解：the radius of the temporal neighborhood used to select points for the temporal exponential filter，即确定平滑的数据，利用哪些数据去做指数平滑，Temporal Smoothing Range：手册的说法是corresponding to the range of the temporal exponential function.但是我没有理解，我觉得这个参数指的是下面这个，可以利用这个参数wise去计算出指数平滑所需要的alpha。具体可以参加这个matlab介绍链接。

3.6 Dialog Box 5 (Space/Time Covariance Analysis)

假设Covariance只和time lag以及spatial lag有关，目标是要拟合出covariance model c(r,t)，要想拟合出c(r,t)，第一步是要生成数据点，也就是不同的time lag和spatial lag下的covariance，生成方式详见手册。第二步就是要拟合Covariance Model，公式如下图所示：

唯一需要说明的是$sum(c_{0i})$需要等于或近似等于数据方差，具体原因可以看手册，也很容易理解了，其中$c_{0i}$就是第i个structure的still。“spatial range” and the ‘temporal range”看具体的方差拟合公式就可以知道，不再详细叙述。具体操作可以多拟合几次，看图效果不错了，就可以停止了。

3.7 Dialog Box 6(BME Estimation)

这个部分就是BME估计部分了，可以进行时间维度的估计和空间维度的估计两部分。 空间分布部分得到的特定时间的BME估计的均值和方差地图，时间分布部分得到的是特定观测点的BME估计得到的关于时间的均值和方差的函数。
Maximum Spatial Distance和Maximum Temporal Distance限定BME对某个点进行估计时的时间和空间范围。
Space/Time Metric: A parameter that is used to calculate the space/time distance. The space/Time distance is obtained as (Spatial distance) + (Space/Time Metric) * (Temporal distance) 简单地说就是一个关于距离和时间的量纲调整系数。
Max Hard Data Point和Max Soft Data Point：BME估计中使用的最多的软硬数据数量。
Local Mean: 这个就直接贴手册了，我也没完全理解。Order of the polynomial used to model the mean trend (or drift) along the spatial and temporal axes within the neighborhood of the estimation point. The default setting is “Zero”, which will use a mean trend of zero and corresponds to simple kriging. “Constant” will use a constant local drift, which corresponds to ordinary kriging applied locally around the estimation point. “Linear” will use a local drift that varies linearly along the spatial and temporal axes. “Quadratic” will use a polynomial of order 2, etc. Generally Order>1 corresponds to universal kriging applied locally around the estimation point.

Estimation Parameters (Spatial Distribution) 部分也有一些参数，这些参数主要是用来对空间分布可视化起作用的。
Estimation Time：需要展示空间BME估计的时间；
Number of Estimation Points (X) and (Y)：x轴和y轴的估计栅格个数；
Area of Estimation Grid：进行BME估计的边界；
mask文件约束估计范围，使得估计范围更加精准，但是也可以不提供。
Estimation Parameters (Temporal Distribution) 也是类似的。
最后绘出的图是这样的：A plot of the time series is displayed when clicking on the tab corresponding to a specific PlotID. The blue solid line displays the BME mean estimates and the green dotted line shows the lower and upper bounds of the 69% confidence interval (which corresponds to the BME mean estimate ± 1 standard deviation under the assumption of a Gaussian distribution). The blue dots show the hard data, while the red triangles and squares show the hardened soft interval and soft Gaussian data. BMEGUI also displays the shape (i.e. either interval or Gaussian) of the PDF describing the soft datum at each soft data point. “Plot List” displays all the estimated time series plots and each entry on the list shows its “Plot ID” and “Station ID”.

至此软件也介绍完了，之后的话可以跑一下软件附带的数据集例子感受一下，也可以熟悉一下。有时间的话建议还是要自己看一遍手册，也比较快，大半天应该就可以看完了。

参考文献

[1]张楚天. 贝叶斯最大熵方法时空预测关键问题研究与应用[D]. 华中农业大学, 2016.
[2]杨勇, 张楚天, 贺立源. 基于贝叶斯最大熵的多因子空间属性预测新方法[J]. 浙江大学学报(农业与生命科学版), 2013, 39(6):636-644.
[3]杨勇, 张若兮. 贝叶斯最大熵地统计方法研究与应用进展[J]. 土壤, 2014, 046(003):402-406.
[4] He J , Kolovos A . Bayesian maximum entropy approach and its applications: a review[J]. Stochastic Environmental Research & Risk Assessment, 2017.