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里面去:
中心极限定理的证明的全部材料都在下面了,你看一看就知道了。一、例子
[例1] 高尔顿钉板试验.
图中每一个黑点表示钉在板上的一颗钉子.每排钉子等距排列,下一排的每个钉子恰在上一排两相邻钉子之间.假设有排钉子,从入口中处放入小圆珠.由于钉板斜放,珠子在下落过程中碰到钉子后以的概率滚向左边,也以的概率滚向右边.如果较大,可以看到许多珠子从处滚到钉板底端的格子的情形如图所示,堆成的曲线近似于正态分布.
如果定义:当第次碰到钉子后滚向右边,令当第次碰到钉子后滚向左边,令.则是独立的,且
那么由图形知小珠最后的位置的分布接近正态.可以想象,当越来越大时接近程度越好.由于时,.因此,显然应考虑的是的极限分布.历史上德莫佛第一个证明了二项分布的极限是正态分布.研究极限分布为正态分布的极限定理称为中心极限定理.
二、中心极限定理
设是独立随机变量序列,假设存在,若对于任意的,成立
称服从中心极限定理.
[例2] 设服从中心极限定理,则服从中心极限定理,其中为数列.
解:服从中心极限定理,则表明
其中.由于,因此
故服从中心极限定理.
三、德莫佛-拉普拉斯中心极限定理
在重贝努里试验中,事件在每次试验中出现的概率为为次试验中事件出现的次数,则
[例3] 用频率估计概率时的误差估计.
由德莫佛—拉普拉斯极限定理,
由此即得
第一类问题是已知,求,这只需查表即可.
第二类问题是已知,要使不小于某定值,应至少做多少次试验?这时利用求出最小的.
第三类问题是已知,求.
解法如下:先找,使得.那么,即.若未知,则利用,可得如下估计: .
[例4] 抛掷一枚均匀的骰子,为了至少有0.95的把握使出现六点的概率与之差不超过0.01,问需要抛掷多少次?
解:由例4中的第二类问题的结论,.即.查表得.将代入,便得. 由此可见,利用比利用契比晓夫不等式要准确得多.
[例5] 已知在重贝努里试验中,事件在每次试验中出现的概率为为次试验中事件出现的次数,则服从二项分布:
的随机变量.求.
解:
因为很大,于是
所以
利用标准正态分布表,就可以求出的值.
[例6] 某单位内部有260架电话分机,每个分机有0.04的时间要用外线通话,可以认为各个电话分机用不用外线是是相互独立的,问总机要备有多少条外线才能以0.95的把握保证各个分机在使用外线时不必等候.
解:以表示第个分机用不用外线,若使用,则令否则令.则.
如果260架电话分机同时要求使用外线的分机数为,显然有.由题意得,
查表得,,故取.于是
取最接近的整数,所以总机至少有16条外线,才能有0.95以上的把握保证各个分机在使用外线时不必等候.
[例7] 根据孟德尔遗传理论,红黄两种番茄杂交第二代结红果植株和结黄果植株的比率为3:1,现在种植杂交种400株,试求结黄果植株介于83和117之间的概率.
解:将观察一株杂交种的果实颜色看作是一次试验,并假定各次试验是独立的.在400株杂交种中结黄果的株数记为,则.
由德莫佛—拉普拉斯极限定理,有
其中,即有
四、林德贝格-勒维中心极限定理
若是独立同分布的随机变量序列,假设,则有
证明:设的特征函数为,则
的特征函数为
又因为,所以
于是特征函数的展开式
从而对任意固定的,有
而是分布的特征函数.因此,
成立.
[例8] 在数值计算时,数用一定位的小数来近似,误差.设是用四舍五入法得到的小数点后五位的数,这时相应的误差可以看作是上的均匀分布.
设有个数,它们的近似数分别是,.,.令
用代替的误差总和.由林德贝格——勒维定理,
以,上式右端为0.997,即以0.997的概率有
[例9] 设为独立同分布的随机变量序列,且互相独立,其中,证明:的分布函数弱收敛于.
证明:为独立同分布的随机变量序列,且互相独立,所以仍是独立同分布的随机变量序列,易知有
由林德贝格——勒维中心极限定理,知的分布函数弱收敛于,结论得证.
作业:
P222 EX 32,33,34,35
五、林德贝尔格条件
设为独立随机变量序列,又
令,对于标准化了的独立随机变量和
的分布
当时,是否会收敛于分布?
[例10] 除以外,其余的均恒等于零,于是.这时就是的分布函数.如果不是正态分布,那么取极限后,分布的极限也就不会是正态分布了.因而,为了使得成立,还应该对随机变量序列加上一些条件.从例题中看出,除以外,其余的均恒等于零,在和式中,只有一项是起突出作用.由此认为,在一般情形下,要使得收敛于分布,在的所有加项中不应该有这种起突出作用的加项.因为考虑加项个数的情况,也就意味着它们都要“均匀地小”.
设是独立随机变量序列,又,,这时
(1)若是连续型随机变量,密度函数为,如果对任意的,有
(2)若是离散型随机变量,的分布列为
如果对于任意的,有
则称满足林德贝尔格条件.
[例11] 以连续型情形为例,验证:林德贝尔格条件保证每个加项是“均匀地小”.
证明: 令,则
于是
从而对任意的,若林德贝尔格条件成立,就有
这个关系式表明, 的每一个加项中最大的项大于的概率要小于零,这就意味着所有加项是“均匀地小”.
六、费勒条件
设是独立随机变量序列,又,,称条件为费勒条件.
林德贝尔格证明了林德贝尔格条件是中心极限定理成立的充分条件,但不是必要条件.费勒指出若费勒条件得到满足,则林德贝尔格条件也是中心极限定理成立的必要条件.
七、林德贝尔格-费勒中心极限定理
引理1 对及任意的,
证明:记,设,由于
因此, ,其次,对,
用归纳法即得.
由于,因此,对也成立.
引理2 对于任意满足及的复数,有
证明:显然
因此,
由归纳法可证结论成立.
引理3 若是特征函数,则也是特征函数,特别地
证明 定义随机变量
其中相互独立,均有特征函数,服从参数的普哇松分布,且与诸 独立,不难验证的特征函数为,由特征函数的性质即知 成立.
林德贝尔格-费勒定理
定理 设为独立随机变量序列,又 .令 ,则
(1)
与费勒条件成立的充要条件是林德贝尔格条件成立.
证明:(1)准备部分
记
(2)
显然(3)
(4)
以及分别表示的特征函数与分布函数,表示的分布函数,那么 (5)
这时
因此林德贝尔格条件化为:对任意,
(6)
现在开始证明定理.设是任意固定的实数.
为证(1)式必须证明
(7)
先证明,在费勒条件成立的假定下,(7)与下式是等价的:
(8)
事实上,由(3)知,又因为
故对一切,
把在原点附近展开,得到
因若费勒条件成立,则对任意的,只要充分大,均有
(9)
这时
(10)
对任意的,只要充分小,就可以有
(11)
因此,由引理3,引理2及(10),(11),只要充分大,就有
(12)
因为可以任意小,故左边趋于0,因此,证得(7)与(8)的等价性.
(2)充分性
先证由林德贝尔格条件可以推出费勒条件.事实上,
(13)
右边与无关,而且可选得任意小对选定的,由林德贝尔格条件(6)知道第二式当足够大时,也可以任意地小,这样,费勒条件成立.
其次证明林德贝尔格条件能保证(1)式成立.注意到(3)及(4),可知,
当时,
当时,
因此
(14)
对任给的,由于的任意性,可选得使,对选定的,用林德贝尔格条件知只要充分大,也可使.因此,已证得了(8),但由于已证过费勒条件成立,这时(8)与(7)是等价的,因而(7)也成立.
(3)必要性
由于(1)成立,因此相应的特征函数应满足(7).但在费勒条件成立时,这又推出了(8),因此,
(15)
上述被积函数的实部非负,故
而且
(16)
因为对任意的,可找到,使,这时由(15),(16)可得
故林德贝尔格条件成立.
八、李雅普诺夫定理
设为独立随机变量序列,又.令,若存在,使有
则对于任意的,有