121 人关注0 条评论
写回答
查看全部 5 个回答
写回答
叶山Shan Ye
GIS/地质/人文地理/可持续发展
A2A 谢邀，
我和我认识的一些人，刚开始用R做空间分析的时候，也遇到过这个问题。R这种开源的东西，优点是各种包很丰富，缺点是有些包的说明写得很乱，地理加权回归（GWR）的R包其实功能很强大，但大部分说明都不大靠谱。
GWR在R里面可以用好几个不同的包来实现，其中步骤最简单的是spgwr。思路就两步：建立窗口、用窗口扫全局。这其实就是GWR本质上的两步。比如我要在全美国范围内统计某两个（或多个）变量之间的回归关系，我可以做一个全局回归（global regression），但因为这些变量在空间分布上或许会有异质性（heterogeneity），表现在统计结果上就是空间不稳定性（nonstationarity），因此只看全局的统计，可能看不出什么结果来。举个不完全恰当但是很容易领会精神的例子，你比如说，我要分析亚洲范围内，经济发展程度与牛肉销量之间的关系，经济越发达的地方，人们就越吃得起牛肉。可是等我统计到印度的时候，坏了，印度大部分人不吃牛肉，这不是经济状况导致的，这一下就影响了全局统计的参考价值，那怎么办呢？我们可以建立一个窗口（正规说法是带宽窗口，bandwidth window），每次只统计窗口范围内的经济与牛肉销量的关系，然后用这个窗口去扫过全局的范围。等统计到印度的时候，印度内部的各地和印度自己比，吃牛肉的人的比例就不会突然减少，这样就能减少这种空间不稳定性对全局统计的影响。
所以，第一步就是要建立这样一个『窗口』。当然了，首先要安装包，我们要用到的R包有：
library(spgwr)
library(rgdal)
library(sf)
library(spData)
library(sp)
library(lattice)
library(ggplot2)
library(ggthemes)
其中，spgwr是做GWR的包，rgdal是用来读取矢量要素的，sf，sp和spData都是用来处理矢量数据的，别的基本都是画图用。
以下默认你会R和GWR的基本操作。并且，以下只展现方法，不要纠结我的数据和结果，我随便找的数据，这个数据本身没有什么意义，所以做出的统计看起来很『壮观』。
我们先导入数据。这里我用的是美国本土48州各个县（county，也有翻译成郡的）的人口普查数据和农业数据，来源是ESRI Online数据库。为啥用这个数据呢？因为...我电脑里面就存了这么个可以用来做GWR的数据...
我们用rgdal读取数据，然后把它画出来看看
require(rgdal)
usa_agri <- readOGR(dsn = "~/Documents/Spatial", layer = "usa_counties")
plot(usa_agri)
会得到这个东西：
readOGR里面，dsn后面加储存shp的路径（加到文件夹为止），layer后面写shp的文件名（不加.shp）。不喜欢rgdal的同学可以不用，用maptools或者spData等别的处理shp的R包代替。不过如果用maptools，要注意处理一下参考系。
我们看一下这个shp里面的列联表都有什么：
可见，shp里面有3108个县的数据，数据有61种。然后再看data下面有什么：
总之就是各种人口普查的数据，后面截不完图，还有经济、房地产和农业之类的数据。那我们就随便选两个来当变量。我就随便挑了，因变量选AVESIZE12，即2012年各个县农场的平均占地面积。自变量选POP_SQMI，也就是人口密度（每平方英里的人口）。
现在正式建立窗口，调用的是spgwr里面的gwr.sel函数：
bw <- gwr.sel( AVE_SIZE12 ~ POP_SQMI, data = usa_agri, gweight = gwr.Gauss,
verbose = FALSE, method = "cv")
其中~前后分别是因变量和自变量。GWR里因变量只能有1个，但自变量可以选多个，如果需要多个自变量的话，就在代码POP_SQMI之后用+号连接就行。gweight是你的空间加权的函数（随空间距离增大而不断衰减的函数，衰减率由下面要提到的带宽控制），这里用的是比较常用的高斯函数，其余的还有gwr.bisquare等函数可以调用。verbose决定是否汇报制定窗口的过程。method是决定构建带宽窗口模型的方法，这里用的cv指的是cross validation，即交叉验证法，也是最常用的方法，简单说就是把数据分成不同的组，分别用不同的方法来做回归计算，计算完了之后记录下结果，然后打乱重新分组，再回归计算，再看结果，周而复始，最后看哪种计算方法的结果最靠谱，这种方法就是最优解。还有一种很常见的选择最佳拟合模型的方法是AIC optimisation法，把method后面的cv改成aic就可以用。具体AIC optimisation是什么：AIC（赤池信息准则）_百度百科。总之，空间加权函数和带宽窗口构建方法的选择是GWR里面十分重要的步骤。
以上便是固定带宽窗口的示意图。比如我在对佐治亚做GWR，这一轮的regression target是红色的这个县，根据做出来的窗口，圆圈以内的县都要被算为红色县的邻县，其权重根据高斯函数等空间权重函数来赋值，而圆圈以外的县，空间权重都赋为0。
不喜欢固定带宽窗口的同学也可以不用它，而是用符合Tobler地理学第一定律的非固定带宽邻域统计，操作方法是在gwr.sel里面加一个命令adapt = TRUE，这样的情况下，根据你设置的k邻居数，每一轮统计的时候，和本轮对象在k以内相邻的多边形的权重参数会被赋值为0到1之间的一个数，比如下图:
我在对佐治亚做GWR，这一轮的regression target是红色的这个县，那么图上标为1的县就是红色县的1阶邻县，标为2的是2阶（邻县的邻县），标为3的是3阶（邻县的邻县的邻县）。如果用非固定带宽邻域统计，k为3，那么1、2、3都被定义为红色县的邻县，它们的权重从3到1依次增加，会按比例被赋上0和1之间的值，而其它没有标注的县，权重为0。
下一步就是用前一步做出的窗口去扫过全局区域：
gwr_result <- gwr(AVE_SIZE12 ~ POP_SQMI, data = usa_agri, bandwidth = bw,
gweight = gwr.Gauss, hatmatrix = TRUE)
这一步如果数据量大，可能会要跑一阵，跑完之后我们看看结果里面有什么：
Call:
gwr(formula = AVE_SIZE12 ~ POP_SQMI, data = usa_agri, bandwidth = bw,
gweight = gwr.Gauss, hatmatrix = TRUE)
Kernel function: gwr.Gauss
Fixed bandwidth: 205880.3
Summary of GWR coefficient estimates at data points:
Min. 1st Qu. Median 3rd Qu.Max. Global
X.Intercept. 7.3883e+01 2.1081e+02 3.2802e+02 6.6691e+02 8.5705e+03 625.5656
POP_SQMI -8.0085e+01 -4.5983e-01 -1.4704e-01 -7.3703e-02 -2.1859e-03 -0.0426
Number of data points: 3108
Effective number of parameters (residual: 2traceS - traceS'S): 119.6193
Effective degrees of freedom (residual: 2traceS - traceS'S): 2988.381
Sigma (residual: 2traceS - traceS'S): 1048.78
Effective number of parameters (model: traceS): 84.90185
Effective degrees of freedom (model: traceS): 3023.098
Sigma (model: traceS): 1042.741
Sigma (ML): 1028.4
AICc (GWR p. 61, eq 2.33p. 96, eq. 4.21): 52109.55
AIC (GWR p. 96, eq. 4.22): 52017.7
Residual sum of squares: 3287040139
Quasi-global R2: 0.4829366
基本上你做GWR该需要的结果这里都有了。比如窗口大小（Fixed bandwidth）是205880.3，意思是前一步构建的带宽窗口是半径205.88千米的圆。Effective number of parameters显示的是你带宽窗口的大小合不合适。Sigma是残差的标准差，这个值要尽量小。Residual sum of squares（RSS）也是对拟合程度的一个评估值。最重要的是最后那个R2，越靠近1说明统计的拟合度越好。我这里面Sigma很大，R2也不是很大，因为我这里只是呈现方法，用的数据本来就是互不相干、没什么太大意义的，所以不用太纠结。如果你是真正的统计数据要来做GWR，就需要注意这些值了。
然后，我们就可以把每个县的R2画在地图上。首先，前面报告里的这些数据，比如R2，要先自己去生成的GWR结果里面去找，然后自己再算一下每个县的local R2，并把它们赋值到shp里面去：
智能运维(AIops)是目前 IT 运维领域最火热的词汇，全称是 Algorithmic IT operations platforms，正规翻译是『基于算法的 IT 运维平台』，直观可见算法是智能运维的核心要素之一。本文主要谈算法对运维的作用，涉及异常检测和归因分析两方面，围绕运维系统Kale 中 skyline、Oculus 模块、Opprentice 系统、Granger causality（格兰杰因果关系）、FastDTW 算法等细节展开。
一、异常检测
异常检测，是运维工程师们最先可能接触的地方了。毕竟监控告警是所有运维工作的基础。设定告警阈值是一项耗时耗力的工作，需要运维人员在充分了解业务的前提下才能进行，还得考虑业务是不是平稳发展状态，否则一两周改动一次，运维工程师绝对是要发疯的。
如果能将这部分工作交给算法来解决，无疑是推翻一座大山。这件事情，机器学习当然可以做到。但是不用机器学习，基于数学统计的算法，同样可以，而且效果也不差。
异常检测之Skyline异常检测模块
2013年，Etsy 开源了一个内部的运维系统，叫 Kale。其中的 skyline 部分，就是做异常检测的模块， 它提供了 9 种异常检测算法 ：
first_hour_average、
simple_stddev_from_moving_average、
stddev_from_moving_average、
mean_subtraction_cumulation、
least_squares
histogram_bins、
grubbs、
median_absolute_deviation、
Kolmogorov-Smirnov_test
简要的概括来说，这9种算法分为两类：
从正态分布入手：假设数据服从高斯分布，可以通过标准差来确定绝大多数数据点的区间；或者根据分布的直方图，落在过少直方里的数据就是异常；或者根据箱体图分析来避免造成长尾影响。
从样本校验入手：采用 Kolmogorov-Smirnov、Shapiro-Wilk、Lilliefor 等非参数校验方法。
这些都是统计学上的算法，而不是机器学习的事情。当然，Etsy 这个 Skyline 项目并不是异常检测的全部。
首先，这里只考虑了一个指标自己的状态，从纵向的时序角度做异常检测。而没有考虑业务的复杂性导致的横向异常。其次，提供了这么多种算法，到底一个指标在哪种算法下判断的更准？这又是一个很难判断的事情。
问题一： 实现上的抉择。同样的样本校验算法，可以用来对比一个指标的当前和历史情况，也可以用来对比多个指标里哪个跟别的指标不一样。
问题二： Skyline 其实自己采用了一种特别朴实和简单的办法来做补充——9 个算法每人一票，投票达到阈值就算数。至于这个阈值，一般算 6 或者 7 这样，即占到大多数即可。
异常检测之Opprentice系统
作为对比，面对相同的问题，百度 SRE 的智能运维是怎么处理的。在去年的 APMcon 上，百度工程师描述 Opprentice 系统的主要思想时，用了这么一张图：
Opprentice 系统的主体流程为：
KPI 数据经过各式 detector 计算得到每个点的诸多 feature；
通过专门的交互工具，由运维人员标记 KPI 数据的异常时间段；
采用随机森林算法做异常分类。
其中 detector 有14种异常检测算法，如下图：
我们可以看到其中很多算法在 Etsy 的 Skyline 里同样存在。不过，为避免给这么多算法调配参数，直接采用的办法是：每个参数的取值范围均等分一下——反正随机森林不要求什么特征工程。如，用 holt-winters 做为一类 detector。holt-winters 有α，β，γ 三个参数，取值范围都是 [0, 1]。那么它就采样为 (0.2, 0.4, 0.6, 0.8)，也就是 4 ** 3 = 64 个可能。那么每个点就此得到 64 个特征值。
异常检测之
Opprentice 系统与 Skyline 很相似
Opprentice 系统整个流程跟 skyline 的思想相似之处在于先通过不同的统计学上的算法来尝试发现异常，然后通过一个多数同意的方式/算法来确定最终的判定结果。
只不过这里百度采用了一个随机森林的算法，来更靠谱一点的投票。而 Etsy 呢？在 skyline 开源几个月后，他们内部又实现了新版本，叫 Thyme。利用了小波分解、傅里叶变换、Mann-whitney 检测等等技术。
另外，社区在 Skyline 上同样做了后续更新，Earthgecko 利用 Tsfresh 模块来提取时序数据的特征值，以此做多时序之间的异常检测。我们可以看到，后续发展的两种 Skyline，依然都没有使用机器学习，而是进一步深度挖掘和调整时序相关的统计学算法。
开源社区除了 Etsy，还有诸多巨头也开源过各式其他的时序异常检测算法库，大多是在 2015 年开始的。列举如下：
Yahoo! 在去年开源的 egads 库。(Java)
Twitter 在去年开源的 anomalydetection 库。(R)
Netflix 在 2015 年开源的 Surus 库。(Pig，基于PCA)
其中 Twitter 这个库还被 port 到 Python 社区，有兴趣的读者也可以试试。
二、归因分析
归因分析是运维工作的下一大块内容，就是收到报警以后的排障。对于简单故障，应对方案一般也很简单，采用 service restart engineering~ 但是在大规模 IT 环境下，通常一个故障会触发或导致大面积的告警发生。如果能从大面积的告警中，找到最紧迫最要紧的那个，肯定能大大的缩短故障恢复时间(MTTR)。
这个故障定位的需求，通常被归类为根因分析（RCA，Root Cause Analysis）。当然，RCA 可不止故障定位一个用途，性能优化的过程通常也是 RCA 的一种。
归因分析之 Oculus 模块
和异常检测一样，做 RCA 同样是可以统计学和机器学习方法并行的~我们还是从统计学的角度开始。依然是 Etsy 的 kale 系统，其中除了做异常检测的 skyline 以外，还有另外一部分，叫 Oculus。而且在 Etsy 重构 kale 2.0 的时候，Oculus 被认为是1.0 最成功的部分，完整保留下来了。
Oculus 的思路，用一句话描述，就是：如果一个监控指标的时间趋势图走势，跟另一个监控指标的趋势图长得比较像，那它们很可能是被同一个根因影响的。那么，如果整体 IT 环境内的时间同步是可靠的，且监控指标的颗粒度比较细的情况下，我们就可能近似的推断：跟一个告警比较像的最早的那个监控指标，应该就是需要重点关注的根因了。
Oculus 截图如下：
这部分使用的 计算方式有两种：
欧式距离，就是不同时序数据，在相同时刻做对比。假如0分0秒，a和b相差1000，0分5秒，也相差1000，依次类推。
FastDTW，则加了一层偏移量，0分0秒的a和0分5秒的b相差1000，0分5秒的a和0分10秒的b也相差1000，依次类推。当然，算法在这个简单假设背后，是有很多降低计算复杂度的具体实现的，这里就不谈了。
唯一可惜的是 Etsy 当初实现 Oculus 是基于 ES 的 0.20 版本，后来该版本一直没有更新。现在停留在这么老版本的 ES 用户应该很少了。除了 Oculus，还有很多其他产品，采用不同的统计学原理，达到类似的效果。
归因分析之 Granger causality
Granger causality（格兰杰因果关系）是一种算法，简单来说它通过比较“已知上一时刻所有信息，这一时刻 X 的概率分布情况”和“已知上一时刻除 Y 以外的所有信息，这一时刻 X 的概率分布情况”，来判断 Y 对 X 是否存在因果关系。
可能有了解过一点机器学习信息的读者会很诧异了：不是说机器只能反应相关性，不能反应因果性的么？需要说明一下，这里的因果，是统计学意义上的因果，不是我们通常哲学意义上的因果。
统计学上的因果定义是：『在宇宙中所有其他事件的发生情况固定不变的条件下，如果一个事件 A 的发生与不发生对于另一个事件 B 的发生的概率有影响，并且这两个事件在时间上有先后顺序（A 前 B 后），那么我们便可以说 A 是 B 的原因。』
归因分析之皮尔逊系数
另一个常用的算法是皮尔逊系数。下图是某 ITOM 软件的实现：
我们可以看到，其主要元素和采用 FastDTW 算法的 Oculus 类似：correlation 表示相关性的评分、lead/lag 表示不同时序数据在时间轴上的偏移量。
皮尔逊系数在 R 语言里可以特别简单的做到。比如我们拿到同时间段的访问量和服务器 CPU 使用率：
然后运行如下命令：
acc_count<-scale(acc$acc_count,center=T,scale=T)
cpu<-scale(acc$cpuload5,center=T,scale=T)
cor.test(acc_count,cpu)
可以看到如下结果输出：
对应的可视化图形如下：
这就说明网站数据访问量和 CPU 存在弱相关，同时从散点图上看两者为非线性关系。因此访问量上升不一定会真正影响 CPU 消耗。
其实 R 语言不太适合嵌入到现有的运维系统中。那这时候使用 Elasticsearch 的工程师就有福了。ES 在大家常用的 metric aggregation、bucket aggregation、pipeline aggregation 之外，还提供了一种 matrix aggregation，目前唯一支持的 matrix_stats 就是采用了皮尔逊系数的计算，接口文档见：
https://www.elastic.co/guide/en/elasticsearch/reference/current/search-aggregations-matrix-stats-aggregation.html
唯一需要注意的就是，要求计算相关性的两个字段必须同时存在于一个 event 里。所以没法直接从现成的 ES 数据中请求不同的 date_histogram，然后计算，需要自己手动整理一遍，转储回 ES 再计算。
饶琛琳，目前就职日志易，有十年运维工作经验。在微博担任系统架构师期间，负责带领11人的SRE团队。著有《网站运维技术与实践》、《ELKstack权威指南》，合译有《Puppet 3 Cookbook》、《Learning Puppet 4》。在众多技术大会上分享过自动化运维与数据分析相关主题。
1、转到VBE界面，菜单工具-引用中看一下有没有丢失的引用，如果有，重新引用一下即可。特别是出错信息提示一些VBA的基本函数未定义，比如LeftDateRight等函数未定义时，90%的情况丢失引用丢失2、Access中执行jetsql语句时使用的很多函数是Access自带的函数，只允许使用在Access界面下，一旦你使用了其他软件做界面，那么很多本来在access+jetsql环境中能够运行的函数将成为错误根源。
3、Access中编写的自定义函数必须由ACCESS环境支持，在其他环境中根本不能使用。 特别注意：许多在VB代码中可以运行的函数并不一定能嵌入到jetsql语句中 以下列出jetsql中的资料供参考： ODBC标量函数 Microsoft
