Saturday, January 30, 2016

Data: International Soil Moisture Network

Summary

The International Soil Moisture Network is an international cooperation to establish and maintain a global in-situ soil moisture database. This database is an essential means of the geoscientific community for validating and improving global satellite observations and land surface models.
Soil moisture, which is the water stored in the upper soil layer, is a crucial parameter for a large number of applications, including numerical weather prediction, flood forecasting, agricultural drought assessment, water resources management, greenhouse gas accounting, civil protection, and epidemiological modeling of water borne diseases. Therefore, the societal benefits of the International Soil Moisture Network are expected to be large.

References

Friday, January 29, 2016

Geography: 地形起伏度

Summary

地形起伏度从数量上讲是单位面积内之高差,或称统计单元内之起伏高度,并得出起伏高度随面积的变化曲线为逻辑斯迪克型。在全国范围内普遍采样的基础上分析得出,中国存在2/10/16/20/22 km2五种不同规模的地形起伏度最佳统计单元,并各具有适用范围。
地形起伏度公式如下
\[LE{R_i} = {E_{\max }} - {E_{\min }}\]
式中:LERi表示以第i个栅格为中心的窗口内的相对高差值,EmaxEmin分别表示该窗口内的最大、最小高程值。窗口分析法的基本原理是:对栅格数据系统中的一个栅格开辟具有固定分析半径的窗口,并在窗口内进行一系列统计计算,窗口在栅格数据矩阵中连续移动完成整个区域的计算。
通常栅格数据的统计单元以矩形为最常见,并且往往是正方形,此类情况又分为奇数和偶数边长。下面分别介绍两种情况下统计单元的计算过程:
Fig. 1
均值变点分析法
曲线上的平均地势起伏度由陡变缓的一点所对应的面积既是最佳统计单元的面积。
将实验区渐变窗口下的平均起伏度值作为均值变点分析法的非线性系统的输出数据,即{Xt, t=1,2,…,n}。令i=2, …, n,对每个i将样本分为两段:X1, X2, …, Xi-1Xi, Xi+1, …, Xn,计算每段样本的算术平均值Xi1Xi2及统计量Si
\[{S_i} = \sum\limits_{t = 1}^{i - 1} {{{\left( {{X_t} - {X_{i1}}} \right)}^2}} + \sum\limits_{t = i}^N {{{\left( {{X_t} - {X_{i2}}} \right)}^2}} \] \[X = \frac{{\sum\limits_{t = 1}^N {{X_t}} }}{N}\] \[S = \sum\limits_{t = i}^N {{{\left( {{X_t} - X} \right)}^2}} \]
最佳统计面积即为S-Si差值最大者所对应的统计面积。

Example

References

[1] 张伟, 李爱农. 基于DEM 的中国地形起伏度适宜计算尺度研究. 地理与地理信息科学. 2012, 28(4): 8~12.
[2] 韩海辉, 高婷, 易欢, 等. 基于变点分析法提取地势起伏度——以青藏高原为例. 地理科学, 2012, 32(1): 101~104.
[3] 张军, 李晓东, 陈春艳, 等. 新疆地势起伏度的分析研究. 兰州大学学报(自然科学版). 2008, 44: 10~13.
[5] 常直杨, 王建, 白世彪, 等. 均值变点分析法在最佳集水面积阈值确定中的应用. 南京师大学报(自然科学版), 2014, 37(1): 147~150.

Wednesday, January 27, 2016

Meteorology: 土壤干燥度

Summary

土壤干燥度的计算采用修正的谢良尼诺夫模型,具体可表为:
\[D = 0.16 \cdot \frac{{{T_{ \ge 10^\circ {\rm{C}}}}}}{{{P_{ \ge 10^\circ {\rm{C}}}}}}\]
式中:D——干燥度;T——年大于10℃积温;P——全年大于10℃降水量之和,在干旱地区也可以近似约等于全年降水量。

Example

References

[1] 胡云锋, 阿拉腾图雅, 艳燕, 于国茂. 内蒙古锡林郭勒生态系统综合监测与评估. 北京: 中国环境出版社, 2013. pp:213.
[2] 张国平, 张增祥, 刘纪远. 中国土壤风力侵蚀空间格局及驱动因子分析. 地理学报, 2001, 56(2):146~158.
[3] 曹祥会, 雷秋良, 龙怀玉, 等. 河北省土壤温度与干湿状况的时空变化特征. 土壤学报, 2015, 52(3):528~537.

Meteorology: 风场强度

Summary

风场强度是影响影响土壤颗粒能否被风力搬运的决定性因素。风场强度的大小可以利用对常规地面气象观测中有关风速、风向的观测数据进一步计算得到。
风场强度的计算采用美国农业部土壤风力侵蚀方程RWEQ中的相应公式:
\[W = \frac{{{n_d}}}{{500}} \cdot \sum\limits_{i = 1}^n {U \cdot {{\left( {U - {U_c}} \right)}^2}} \]
式中:W——风能强度因子,m3/s3n——观测次数;nd——观测活动包含天数;U——离地面2 m高处的风速;Uc——2 m高处的临界风速,一般设置为5 m/s。

Example

Practice.rar.包含ANUSPLIN批量并行插值代码

References

[1] 胡云锋, 阿拉腾图雅, 艳燕, 于国茂. 内蒙古锡林郭勒生态系统综合监测与评估. 北京: 中国环境出版社, 2013. pp:211~213.
[2] 巩国丽, 刘纪远, 邵全琴. 基于RWEQ 的20 世纪90 年代以来内蒙古锡林郭勒盟土壤风蚀研究. 地理科学进展, 2014, 33(6): 825~834.
[3] Feras Youssef, Saskia Visser, Derek Karssenberg, et al.. Calibration of RWEQ in a patchy landscape; a first step towards a regional scale wind erosion model. Aeolian Research, 2012, (3): 467~476.

Saturday, January 23, 2016

Matlab: 数据类型转换

Summary

数据类型的转换需要考虑原类型及目标类型对应数据的范围,比如:原数值包含负数,在不改变数值的前提下则无符号数据类型则不可以使用了。当然,实际情况往往多种多样,具体问题还需具体分析。

Example

示例数据的数值类型是single,并且包含负值(-99),将此转换为整形,考虑到数值的范围,因此int16比较合适。
代码运行流畅,结果正确。

References

Thursday, January 21, 2016

Math: 灰靶决策

Fig. 1
靶心的设置直接决定着评价结果。如Fig. 1所示,A~C中心的红点即为靶心,靶心指最优或最适状态,此处它是无大小且不可再分的图形,满足数学定义的点,A表示最简单、最简洁的灰靶,评价靶心仅为一个数值,评价的过程就是量测指标与靶心之间的距离,这类指标如NPP、植被覆盖度,进一步,在A的基础上得到的评价结果还可以划分为若干等级,绿色、黄色和蓝色就对应不同的类别,还有一类指标不能仅依靠唯一数值作为其评价的靶心,它的靶心应设置为一个区间,如C在红点之外围还有一圈红色区域围绕,此时的靶心不再是唯一的数值,而是一个范围,指标在此范围就是位于靶心位置。File.cdr

Friday, January 15, 2016

Math: 可拓综合评价法的建模步骤

矛盾问题,就是指在现有条件下无法实现人们要达到的目标的问题。

Summary

利用可拓综合评价法,一般要经历确定经典域、确定节域、确定待评价物元、确定评价指标的权重、确定待评价事物关于各类别的关联度以及确定待评价事物的类别和级别变量特征值六个步骤。

Step 1 确定经典域

设有m个评价物元(或评价类别)N1N2,…,Nm,将各评价物元(或类别)对应的特征值范围用[aij, bij]表示,则同征物元R0可表示为:
\[{R_0} = \left[ {\begin{array}{*{20}{c}} N&{{N_1}}&{{N_2}}& \cdots &{{N_m}}\\ {{c_1}}&{\left[ {{a_{11}},{b_{11}}} \right]}&{\left[ {{a_{12}},{b_{12}}} \right]}& \cdots &{\left[ {{a_{1m}},{b_{1m}}} \right]}\\ {{c_2}}&{\left[ {{a_{21}},{b_{21}}} \right]}&{\left[ {{a_{22}},{b_{22}}} \right]}& \cdots &{\left[ {{a_{2m}},{b_{2m}}} \right]}\\ \vdots & \vdots & \vdots & \vdots & \vdots \\ {{c_n}}&{\left[ {{a_{n1}},{b_{n1}}} \right]}&{\left[ {{a_{n2}},{b_{n2}}} \right]}& \cdots &{\left[ {{a_{nm}},{b_{nm}}} \right]} \end{array}} \right]\]
式中:Nj——第j个评价物元(或类别);ci——第i个评价指标;vij=[aij,bij]——Nj关于ci所规定的量值范围,即经典域。

Step 2 确定节域

\[{R_p} = \left( {P,C,{V_p}} \right) = \left[ {\begin{array}{*{20}{c}} P&{{c_1}}&{\left[ {{a_{1p}},{b_{1p}}} \right]}\\ {}&{{c_2}}&{\left[ {{a_{2p}},{b_{2p}}} \right]}\\ {}& \vdots & \vdots \\ {}&{{c_n}}&{\left[ {{a_{np}},{b_{np}}} \right]} \end{array}} \right]\]
式中:P——评价类别的全体;[aipbip]——P关于ci所取的量值范围,即节域。

Step 3 确定待评价物元

对待评价物元,把监测数据或分析结果用物元表示为:
\[{R_d} = \left[ {\begin{array}{*{20}{c}} P&{{c_1}}&{{v_1}}\\ {}&{{c_2}}&{{v_2}}\\ {}& \vdots & \vdots \\ {}&{{c_n}}&{{v_n}} \end{array}} \right]\]
式中:Rd——待评价物元;vi——待评价事物对应于ci的数值。

Step 4 确定评价指标的权重

\[{W_i} \ge 0\left( {i = 1,2, \cdots ,n} \right)\] \[\sum\limits_{i = 1}^n {{W_i}} = 1\]

Step 5 确定关联度

计算距:
\[\begin{array}{c} \rho \left( {{v_i},{v_{ij}}} \right) = \left| {{v_i} - \frac{1}{2}\left( {{a_{ij}} + {b_{ij}}} \right)} \right| - \frac{1}{2}\left( {{b_{ij}} - {a_{ij}}} \right)\\ \rho \left( {{v_i},{v_{ip}}} \right) = \left| {{v_i} - \frac{1}{2}\left( {{a_{ip}} + {b_{ip}}} \right)} \right| - \frac{1}{2}\left( {{b_{ip}} - {a_{ip}}} \right) \end{array}\]
式中:ρvi,vij)——点vi与区间vij的距;ρvi,vip)——点vi与区间vip的距。
计算关联函数:
\[{K_j}\left( {{v_i}} \right) = \left\{ {\begin{array}{*{20}{c}} {\frac{{\rho \left( {{v_i},{v_{ij}}} \right)}}{{\rho \left( {{v_i},{v_{ip}}} \right) - \rho \left( {{v_i},{v_{ij}}} \right)}}{v_i} \notin {v_{ij}}}\\ {\frac{{ - \rho \left( {{v_i},{v_{ij}}} \right)}}{{\left| {{v_{ij}}} \right|}}{v_i} \in {v_{ij}}} \end{array}} \right.\]
Kjvi)——关联函数,待评价事物的指标ci关于类别j的归属度;|vij|——区间[aijbji]的长度,即|bij-aij|。
计算关联度:
\[{K_j}\left( p \right) = \sum\limits_{i = 1}^n {{W_i}{K_j}\left( {{v_i}} \right)} \]
式中:Kjp)——在考虑指标权重下,待评价事物各指标ci关于类别j的关联度的组合值。

Step 6 确定待评价事物的类别和级别变量特征值

Kj0p)=maxKjp)(j=1,2,…,m),则评定p属于类别j0。记:
\[\overline {{K_j}} \left( p \right) = \frac{{{K_j}\left( p \right) - \mathop {\min }\limits_{1 \le j \le m} {K_j}\left( p \right)}}{{\mathop {\max }\limits_{1 \le j \le m} {K_j}\left( p \right) - \mathop {\min }\limits_{1 \le j \le m} {K_j}\left( p \right)}}\]
p的级别变量特征值j*为:
\[{j^ * } = \frac{{\sum\limits_{j = 1}^m {j \cdot \overline {{K_j}} \left( p \right)} }}{{\sum\limits_{j = 1}^m {\overline {{K_j}} \left( p \right)} }}\]

References

[1] 何逢标. 综合评价方法MATLAB实现. 北京: 中国社会科学出版社. pp: 269~270.
[2] 方统中, 杜耘, 蔡述明, 等. 模糊数学在洪湖营养化评价中的应用. 浙江林学院学报, 2008, 25(4):517~521.
[3] 高军省, 吕小凡. 湖泊富营养化评价的可拓学方法及应用. 长江大学学报(自然科学版). 2010, 7(1):33~37.

Thursday, January 14, 2016

Math: 可拓综合评价法概述

Summary

模糊综合评价法是建立在模糊集合基础上的评价方法,利用它可解决“内涵明确,外延不明确”的具有模糊性的一类评价问题。而可拓综合评价法是建立在可拓集合基础上的评价方法,它不仅可以从数量上刻画被评价对象本身在状态的所属程度,而且还可以从数量上刻画何时为一种性态与另一种性态的边界。
在日常生活中,我们常说“量变引起质变”,但究竟多少“量”的累积才会引起“质”的变化,通常难以说明。“可拓学是一种把事物的质与量有机结合起来的理论,它从定性与定量两个角度去研究解决具有矛盾的问题,其理论支柱是物元理论和可拓集合论,其逻辑细胞是物元”。
可拓学是我国学者蔡文于1983年创立的一门新学科,它用形式化的模型研究事物拓展的可能性和开拓创新的规律与方法。物元以有序的三元组R=(NCV)来表达。若事物Nn个特征Cc1c2,…,cn),待评价的m个物元R1=(N1CV1),R2=(N2CV2),…,Rm=(NmCVm)具有相同的特征C,则具有相同特征的物元体R可表示为:
\[R = \left[ {\begin{array}{*{20}{c}} N&{{N_1}}&{{N_2}}& \cdots &{{N_m}}\\ {{c_1}}&{{v_{11}}}&{{v_{12}}}& \cdots &{{v_{1m}}}\\ {{c_2}}&{{v_{21}}}&{{v_{22}}}& \cdots &{{v_{2m}}}\\ \vdots & \vdots & \vdots & \vdots & \vdots \\ {{c_n}}&{{v_{n1}}}&{{v_{n2}}}& \cdots &{{v_{nm}}} \end{array}} \right]\]
其中,Ni——待评价的物元;N——待评价物元N1N2,…,Nm的全体;vij——第j个待评价物元第i个特征的量值,i=1,2,…,nj=1,2,…,m
以此为基础,建立的可拓综合评价法,通过具有相同特征的物元的表示、各特征量值范围的确定、待评价物元各项特征值的描述、关联度的计算,最终得出评价物元所属的类别,现今在水利、环境、经济与管理等众多领域得到了越来越多的应用。

References

[1] 何逢标. 综合评价方法MATLAB实现. 北京: 中国社会科学出版社. pp: 269~270.

Tuesday, January 12, 2016

Math: 熵权法

Description

物理学上的熵(entropy)用来表示任何一种能量在空间中的分布的均匀程度,能量分布的越均匀,熵就越大。当一个体系的能量完全均匀分布时,这个系统的熵就达到最大值。在信息论中,熵又称为平均信息量,它是信息的一个度量。
利用熵的概念确定指标权重的方法称为熵权法。其出发点是根据某同一指标观测值之间的差异程度来反映其重要程度,如果各被评价对象的某项指标的数据差异不大,则反映该指标对评价系统所起的作用不大。
熵权法是一种客观的赋权方法,它是利用各指标的熵值所提供的信息量的大小来决定指标权重的方法。熵权法的作用有:用熵权法给指标赋权可以避免各评价指标权重的人为因素干扰,使评价结果更符合实际;通过对各指标熵的计算,可以衡量出指标信息量的大小,从而确保所建立的指标能反映绝大部分的原始信息。

Method

(1)形成决策矩阵
设参与评价的对象集为M=(M1M2,…,Mm),指标集为D=(D1D2,…,Dn),评价对象Mi对指标Dj的值记为xiji=1,2,…,mj=1,2,…,n),则形成的决策矩阵X为:
\[X = \left[ {\begin{array}{*{20}{c}} {}&{{D_1}}&{{D_2}}& \cdots &{{D_n}}\\ {{M_1}}&{{x_{11}}}&{{x_{12}}}& \cdots &{{x_{1n}}}\\ {{M_2}}&{{x_{21}}}&{x{}_{22}}& \cdots &{{x_{2n}}}\\ \vdots & \vdots & \vdots & \vdots & \vdots \\ {{M_m}}&{{x_{m1}}}&{{x_{m2}}}& \cdots &{{x_{mn}}} \end{array}} \right]\]
(2)标准化决策矩阵
为了消除各指标量纲不同对方案决策带来的影响,或者处理一些指标值为负的决策问题,对决策矩阵X进行标准化处理,从而形成标准化矩阵V=(vijm×n。根据指标的性质,将指标分为两类。一类是越大越优型指标,也称为效益型指标;另一类是越小越优型指标,也称为成本型指标。
标准化处理时根据指标性质,采用相应的标准化形式:
对于越大越优型指标:
\[{v_{ij}} = \frac{{{x_{ij}} - \min \left( {{x_j}} \right)}}{{\max \left( {{x_j}} \right) - \min \left( {{x_j}} \right)}}\]
对于越小越优型指标:
\[{v_{ij}} = \frac{{\max \left( {{x_j}} \right) - {x_{ij}}}}{{\max \left( {{x_j}} \right) - \min \left( {{x_j}} \right)}}\]
(3)计算第j项指标下,第i个评价对象的特征比重
对于某一个指标jvij的值差异越大,表明该项指标对于被评价对象的作用越大,即该项指标提供给被评价对象的有用信息越多。根据熵的概念,信息的增加意味着熵的减少,熵可以用来度量这种信息量的大小。
记第j项指标下,第i个评价对象的特征比重为pij,则
\[{p_{ij}} = \frac{{{v_{ij}}}}{{\sum\limits_{i = 1}^m {{v_{ij}}} }}\]
(4)计算第j项指标的熵值ej
\[{e_j} = - \frac{1}{{\ln \left( m \right)}}\sum\limits_{i = 1}^m {{p_{ij}}\ln \left( {{p_{ij}}} \right)} \] 当pij=0或者pij=1时,认为pijln(pij)=0。
(5)计算第j项指标的差异性系数dj
观察熵值的计算公式,对于某一项指标Djvij的差异越小,ej越大。当各被评价对象第j项指标值全部相等时,ej=emax=1。根据熵的概念,各被评价对象第j项指标值差异越大,表明该指标反映的信息量越大。因此,定义差异系数dj
\[{d_j} = 1 - {e_j}\]
(6)确定各指标的熵权
\[{w_j} = \frac{{{d_j}}}{{\sum\limits_{k = 1}^n {{d_k}} }}j = 1,2, \cdots ,n\]

References

[1] 郭亚军. 综合评价理论、方法及应用. 北京: 科学出版社, 2007. Pp15~16.
[2] 何鑫, 朱宏泉, 高成凤. 基于熵权法与TOPSIS法的房地产项目投资风险评价. 商业研究, 2009, 03:105~108.

Monday, January 11, 2016

Math: 评价指标的一致化

Description

一般来说,指标x1, x2, …, xm中,可能含有“极大型”指标、“极小型”指标、“居中型”指标和“区间型”指标。对于某些定量指标,如产值、利润等,我们自然期望它们的取值越大越好,这类指标我们称之为极大型指标;而对于诸如成本、能耗等一类指标,我们自然期望它们的取值越小越好,这类指标称之为极小型指标;诸如人的身高、体重等指标,我们既不期望它们的取值越大越好,也不期望期望它们的取值越小越好,而是期望它们的取值越居中越好,我们称这类指标为居中型指标;而区间型指标是期望其取值以落在某个区间内为最佳的指标。
若指标x1, x2, …, xm中既有极大型指标、极小型指标,又有居中型指标或区间型指标,则必须在对各备选方案进行综合评价之前,将评价指标的类型作一致化处理。否则,就无法定性地判定综合评价函数中的yi值是否是取值越大越好,或是取值越小越好,或是取值越居中越好。因此,也就无法根据y值的大小来综合评价各备选方案的优劣。因此,需将指标作类型一致化的处理。
对于极小型指标x,令
\[{x^*} = M - x\]
\[{x^*} = \frac{1}{x}\left( {x > 0} \right)\]
式中,M为指标x的一个允许上界。
对于居中型指标x,令
\[{x^*} = \left\{ {\begin{array}{*{20}{c}} {2\left( {x - m} \right),m \le x \le \frac{{M + m}}{2}}\\ {2\left( {M - x} \right),\frac{{M + m}}{2} \le x \le M} \end{array}} \right.\]
式中,m为指标x的一个允许下界,M为指标x的一个允许上界。
对于区间型指标x,令
\[{x^*} = \left\{ {\begin{array}{*{20}{c}} {1.0 - \frac{{{q_1} - x}}{{\max \left\{ {{q_1} - m,M - {q_2}} \right\}}},x < {q_1}}\\ {1.0,x \in \left[ {{q_1},{q_2}} \right]}\\ {1.0 - \frac{{x - {q_2}}}{{\max \left\{ {{q_1} - m,M - {q_2}} \right\}}},x > {q_2}} \end{array}} \right.\]
式中,[q1,q2]为指标x的最佳稳定区间,Mm分别为x的允许上、下界。
如此,非极大型评价指标x通过上式可以转换为极大型指标了。

References

[1] 郭亚军. 综合评价理论、方法及应用. 北京: 科学出版社, 2007. Pp15~16.

Geography: 典型草原的自然保护区

Summary

依据ProtectedPlanet网站,典型草原研究区已知有两个自然保护区,分别如图所示范围,锡林郭勒草原国家级自然保护区和乌拉盖湿地自然保护区,分别建立于1987年和2001年。
Fig. 1

References

Wednesday, January 6, 2016

Data: Red List Spatial Data

Summary

The IUCN Red List of Threatened Species contains assessments for just over 76,000 species, of which about two-thirds have spatial data. This spatial data provided below is for comprehensively assessed taxonomic groups. It is important to note that some species such as those listed as Data Deficient are not mapped and subspecies are mapped within the parental species. The data is available as ESRI shapefiles format and contains the known range of each species. Ranges are depicted as polygons, except for the freshwater HydroSHED tables. The shapefiles contain taxonomic information, distribution status, IUCN Red List category, sources and other details about the maps.
Fig. 1

References

Tuesday, January 5, 2016

Data: UN Population Division

Summary

The 2015 Revision of World Population Prospects is the twenty-fourth round of official United Nations population estimates and projections that have been prepared by the Population Division of the Department of Economic and Social Affairs of the United Nations Secretariat. The main results are presented in a series of Excel files displaying key demographic indicators for each development group, income group, major area, region and country for selected periods or dates within 1950-2100. A publication labelled Key findings and advance tables, which provide insights on the results of this latest revision, is also made available here.
Fig. 1

References

Data: ProtectedPlanet

Summary

ProtectedPlanet.net is the online interface for the World Database on Protected Areas (WDPA), a join project of IUCN and UNEP, and the most comprehensive global database on terrestrial and marine protected areas. ProtectedPlanet.net lets you discover the protected areas of the world through exploring the maps and intuitive searching, feeding you information from the WDPA, photos from Panoramio and text descriptions from Wikipedia.
Fig. 1

References

Data: WorldPop

Summary

High spatial resolution, contemporary data on human population distributions are a prerequisite for the accurate measurement of the impacts of population growth, for monitoring changes and for planning interventions. The WorldPop project aims to meet these needs through the provision of detailed and open access population distribution datasets built using transparent approaches.
The WorldPop project was initiated in October 2013 to combine the AfriPop, AsiaPop and AmeriPop population mapping projects. It aims to provide an open access archive of spatial demographic datasets for Central and South America, Africa and Asia to support development, disaster response and health applications. The methods used are designed with full open access and operational application in mind, using transparent, fully documented and peer-reviewed methods to produce easily updatable maps with accompanying metadata and measures of uncertainty.

China

Fig. 1

References

Saturday, January 2, 2016

Tool: Data Structure Visualizations

Summary

理解复杂算法的最好方法就是观察他们计算的过程。这里就有一个查看多种数据结构及算法的交互动画网站。

Example

点击左上方的Algorithms进入算法可视化界面:
Fig. 1
以Sorting(排序)为例,单击Sorting-Comparison Sorting,如Fig. 2,上方红框是各类排序方法,如Bubble Sort冒泡排序等等,下方红框是General Animation Controls动画过程控制按钮。
Fig. 2

References

Friday, January 1, 2016

Matlab+Python: Theil–Sen estimator

Summary

In non-parametric statistics, there is a method for robust simple linear regression that chooses the median slope among all lines through pairs of two-dimensional sample points. It has been called the Theil–Sen estimator, Sen's slope estimator, slope selection, the single median method, the Kendall robust line-fit method, and the Kendall–Theil robust line. It is named after Henri Theil and Pranab K. Sen, who published papers on this method in 1950 and 1968 respectively.

Example

本例以Anscombe's Quartet为例,检验Theil–Sen estimator对趋势的估计效果。关于Anscombe's Quartet四组数据的特征请参考这里的描述。这里的代码整合Matlab+Python,Theil–Sen estimator计算调用Python Scipy。
Fig. 1
观察Fig. 1,蓝线是简单线性回归,红线是Theil–Sen estimator,两条灰虚线之间表示α=0.9水平下的置信区间。A与D,蓝线与红线基本重合,A的蓝线基本可以反映数据的趋势,但D的蓝线显然不能反映数据分布的趋势,此处的红线也还是不能很好的剔除异常值的影响。B组数据并非呈现线性趋势,所以红线并不会比蓝线更好。C组红线非常好地剔除一个异常值的影响,很好地诠释其他10个数据分布趋势。综上,Theil–Sen estimator的确比简单线性回归优秀,但是也同样存在着无效的情况。
以下Fig. 2~5是不同α水平下的置信区间范围,很显然,随着α增大区间范围增大,更多的数据进入置信区间。
Fig. 2
Fig. 3
Fig. 4
Fig. 5
Fig. 6

References