白话空间统计番外篇：中位数中心算法[转]
昨天在介绍中位数中心算法的时候，挖了个巨大的坑，结果导致老夫一夜没有睡好，脑子里面飞来飞去全部都是各种选择和迭代算法，今天终于下定决心把这个坑给填上。
其实我一直是不愿意填算法坑的……主要是自己的数学水平很一般，很容易出现填坑不成自己反被埋的情况，但是这个坑不填又不行，所以在填坑之前说明：这个仅是虾神我自己的理解，不代表原文（限于能力问题，数学论文确实不怎么能读透），如果有疑惑或者错误，请自行查阅原始论文，虾神只负责科普。
地理上有若干个点，现在需要寻找一个到所有点距离总和最小的点，是几何学和区域学研究中的一个非常重要和著名的问题，但是最早由一个经济学家阿尔弗雷德·韦伯（德语：Alfred Weber，1868年6月30日－1958年5月2日，如下图）提出，所以也称之为“韦伯问题（The Weberproblem）”。这位韦伯先生并非是一个数学家，而是一位经济学家、社会学家和文化理论家。阿尔弗雷德·韦伯是马克斯·韦伯（组织理论之父）的弟弟。他创立了工业区位理论，深刻影响了现代经济地理学的发展。
在讲韦伯问题之前，要涉及到数学里面的一个重大的问题，就是所谓的费马点（Fermat Point)的问题。费马点是17世纪的一个法国律师皮耶·德·费马（Pierrede Fermat，如下图）提出来的，这位专职律师被称为“业余数学之王”，在数学的神殿里面，有他一张王座镇压诸天……费马最后定理在中国习惯称为费马大定理，西方数学界原名“最后”的意思是：其它猜想都证实了，这是最后一个。（PS:貌似不务正业的人都很厉害，国外有卖缝纫机的托马斯.沃森（IBM的创始人），国内有教英语的马云……）
费马点是啥东西呢？
费马点，指的是，在三角形内部，有一个点，这个点到三角形三个顶点的距离之和最短，如下图：
如果三角形的三个内角都小于120度，如上图：三角形ABC内部的一个点D，就是离这个三角形三个顶点距离总和最近的一个点，从这个点向三角形三个顶点连线，得出的三个角正好整分费马点所在的周角，即均为120度。所以费马点也称为三角形的等角中心。
那么如果三角形有一个内角大于等于120度，那么这个钝角的顶点就是费马点。
要找到这个费马点，也非常的容易，根本不需要我们去迭代计算（数学是一种追求完美的学科），只需要将三角形的三条边都做一个等边三角形，然后用这个等边三角形做一个外接圆，三个外接圆的交点，正好就是这个费马点。（等边三角形做外接圆的方法就不详说了，太简单了）。
四边形的费马点就更容易了，凸四边形，费马点就是对角线交点；凹四边形，费马点就是凹点。
我们知道，三角形和四边形这种在数学上非常特殊的情况，在现实生活在确实不多，特别是在多点之间计算费马点的话。
所以复杂多边形内的费马点计算，也一直都是数学界津津乐道的话题。
因为复杂多边形内的费马点没有公式来实现，所以到现在为止想通过一个公式就计算出费马点是不可能的事情。当然，初等数论里面，对正多边形提出了一些方法，例如分割成多个三角形这种方法，但是仅适用于正多边形。
而且费马点的寻找是不涉及到任何权重的，所以算出来的结果是完全几何结果，几何图形的费马点完全是在图形内部。
韦伯问题在费马点的基础上扩展了权重概念，那么带来的一个问题就是：有些点的权重，被设置为了负数，结果就是加权计算的时候，会让这个中位数中心可能跳出几何范围之外。
例如：还是仓库运输的问题，我们要计算中位数中心，那么如果有一个仓库不但可以存放中转货物，还提供了加油、修车、保养、司机休息的服务……那么这个仓库的权重计算可能就会被设置为负数，哪怕他的距离可能很大，但是所有来到这个仓库的车辆，都会直接忽略距离因素而达到更优化的效果。
1962年，普林斯顿大学的哈罗德.威廉.库恩和罗伯特.E.库伦提出了迭代最小二乘法解决韦伯问题的方案。具体的算法如下：
首先确定起算点，起算点的确定非常简单，Kuhn和Kuenne选择了最简单的一种起算方法，就是几何平均数，用所有点的平均中心作为起算点。
接下去就是确定迭代优化方案了。迭代的优化方案主要就是对候选点的选择，有如下几个关键：
1、以起算点为参照，在什么地方选择候选点（方向）。
2、候选点选在起算点多远的地方比较好（距离）。
先讲方向的问题，理论上来说，只需要向任意一个不同的方向移动，就可以了。随便向任意方向移动，都会产生不同的距离总和。
生成新的距离总和之后，与原来的起算点的距离总和进行对比，如果大于原来的距离总和，就说明这个候选点是错误的，丢弃，重新寻找。如果小于原来的距离总和，说明比起算点要优化，将他设为新的起算点，也就是候选点，然后以这个新的起算点，迭代进行选择寻找。
然后再讲讲移动多少距离合适。
Kuhn和Kuenne在他们的论文里面，设定了一种很实用的距离公式，就是所谓的Weiszfeld 算法，这是一种重复加权最小二乘法，如下所示：
第一次选择候选点的距离的时候，直接采用所有的点的平均距离y作为移动距离，移动完成之后计算，并且把这个y带入到公式中，求解出下一次需要移动的距离。
根据我们迭代的次数的增加，会发现数据会逐渐的收敛。最后可以计算出最优的候选点，作为最后的位置。
当然，如果有权重的话，在每个点上面，还需要加上权重进行计算，如下公式：
其中Wi就是每个点的权重。
理论上，这个优化可以无限的接近无穷大（与最优点的距离无限接近于0），但是无论是计算结果还是计算机能表示的结果，都是有一个极限的，在实际的计算中，这个极限就是我们可以接受的精度。
一般来说，在GIS里面，就是你创建坐标系和要素图层的时候指定的精度，就是默认接受的精度值。
当然，里面还有很多很多其他的东西，比如各种条件什么的，我这里就不一一说明了，有兴趣的同学，请参考如下文章：
https://en.wikipedia.org/wiki/Geometric_median
