如何找到数据点集群的中心？
Rya*_*yan 32 algorithm geocoding cluster-analysis data-mining markerclusterer
假设我在过去的一年中每天都绘制了一架直升机的位置,并提出了以下地图:
任何看到这个的人都可以告诉我这架直升机是在芝加哥以外的.
如何在代码中找到相同的结果？
我正在寻找这样的东西:
$geoCodeArray = array([GET=http://pastebin.com/grVsbgL9]);
function findHome($geoCodeArray) {
// magic
return $geoCode;
}
最终生成这样的东西:
更新:示例数据集
这是一张包含样本数据集的地图:http://batchgeo.com/map/c3676fe29985f00e1605cd4f86920179
这是一个包含150个地理编码的pastebin:http://pastebin.com/grVsbgL9
以上包含150个地理编码.前50个在靠近芝加哥的几个集群中.其余的分布在全国各地,包括纽约,洛杉矶和旧金山的一些小集群.
我有大约一百万(严重)这样的数据集,我需要迭代并确定最可能的"家".非常感谢您的帮助.
更新2:飞机切换到直升机
飞机概念引起了对物理机场的过多关注.坐标可以在世界的任何地方,而不仅仅是机场.让我们假设它是一架超级直升机,不受物理,燃料或其他任何东西的束缚.它可以降落在它想要的地方.;)
Wal*_*oss 15
通过将纬度和经度转换为笛卡尔坐标,即使点遍布地球,以下解法仍然有效.它执行一种KDE(内核密度估计),但在第一次传递中,仅在数据点处评估内核的总和.应该选择内核以适应问题.在下面的代码中,我可以开玩笑地/擅自称为Trossian,即d≤h为2-d²/h²,d> h为h²/d²(其中d为欧几里德距离,h为"带宽" $global_kernel_radius) ,但也可能是高斯(e -d²/ 2h²),Epanechnikov内核(d <h为1-d²/h²,否则为0)或其他内核.可选的第二遍通过对本地网格上的独立内核求和,或通过计算质心,在本地优化搜索,$local_grid_radius.
本质上,每个点对它周围的所有点(包括其自身)求和,如果它们更接近(通过钟形曲线),则称它们更多,并且还通过可选的重量阵列称重它们$w_arr.获胜者是具有最大金额的点.一旦找到胜利者,我们正在寻找的"家"可以通过在获胜者周围重复相同的过程(使用另一个钟形曲线),或者可以估计为所有点的"质心"在获胜者的给定半径范围内,半径可以为零.
通过选择适当的内核,选择如何在本地优化搜索以及调整参数,算法必须适应问题.对于示例数据集,第一遍的Trossian内核和第二遍的Epanechnikov内核,所有3个半径设置为30英里,网格步长为1英里可能是一个很好的起点,但只有两个子芝加哥群集应被视为一个大集群.否则必须选择较小的半径.
function find_home($lat_arr, $lng_arr, $global_kernel_radius,
$local_kernel_radius,
$local_grid_radius, // 0 for no 2nd pass
$local_grid_step, // 0 for centroid
$units='mi',
$w_arr=null)
{
// for lat,lng <-> x,y,z see http://en.wikipedia.org/wiki/Geodetic_datum
// for K and h see http://en.wikipedia.org/wiki/Kernel_density_estimation
switch (strtolower($units)) {
/* */case 'nm' :
/*or*/case 'nmi': $m_divisor = 1852;
break;case 'mi': $m_divisor = 1609.344;
break;case 'km': $m_divisor = 1000;
break;case 'm': $m_divisor = 1;
break;default: return false;
}
$a = 6378137 / $m_divisor; // Earth semi-major axis (WGS84)
$e2 = 6.69437999014E-3; // First eccentricity squared (WGS84)
$lat_lng_count = count($lat_arr);
if ( !$w_arr) {
$w_arr = array_fill(0, $lat_lng_count, 1.0);
}
$x_arr = array();
$y_arr = array();
$z_arr = array();
$rad = M_PI / 180;
$one_e2 = 1 - $e2;
for ($i = 0; $i < $lat_lng_count; $i++) {
$lat = $lat_arr[$i];
$lng = $lng_arr[$i];
$sin_lat = sin($lat * $rad);
$sin_lng = sin($lng * $rad);
$cos_lat = cos($lat * $rad);
$cos_lng = cos($lng * $rad);
// height = 0 (!)
$N = $a / sqrt(1 - $e2 * $sin_lat * $sin_lat);
$x_arr[$i] = $N * $cos_lat * $cos_lng;
$y_arr[$i] = $N * $cos_lat * $sin_lng;
$z_arr[$i] = $N * $one_e2 * $sin_lat;
}
$h = $global_kernel_radius;
$h2 = $h * $h;
$max_K_sum = -1;
$max_K_sum_idx = -1;
for ($i = 0; $i < $lat_lng_count; $i++) {
$xi = $x_arr[$i];
$yi = $y_arr[$i];
$zi = $z_arr[$i];
$K_sum = 0;
for ($j = 0; $j < $lat_lng_count; $j++) {
$dx = $xi - $x_arr[$j];
$dy = $yi - $y_arr[$j];
$dz = $zi - $z_arr[$j];
$d2 = $dx * $dx + $dy * $dy + $dz * $dz;
$K_sum += $w_arr[$j] * ($d2 <= $h2 ? (2 - $d2 / $h2) : $h2 / $d2); // Trossian ;-)
// $K_sum += $w_arr[$j] * exp(-0.5 * $d2 / $h2); // Gaussian
}
if ($max_K_sum < $K_sum) {
$max_K_sum = $K_sum;
$max_K_sum_i = $i;
}
}
$winner_x = $x_arr [$max_K_sum_i];
$winner_y = $y_arr [$max_K_sum_i];
$winner_z = $z_arr [$max_K_sum_i];
$winner_lat = $lat_arr[$max_K_sum_i];
$winner_lng = $lng_arr[$max_K_sum_i];
$sin_winner_lat = sin($winner_lat * $rad);
$cos_winner_lat = cos($winner_lat * $rad);
$sin_winner_lng = sin($winner_lng * $rad);
$cos_winner_lng = cos($winner_lng * $rad);
$east_x = -$local_grid_step * $sin_winner_lng;
$east_y = $local_grid_step * $cos_winner_lng;
$east_z = 0;
$north_x = -$local_grid_step * $sin_winner_lat * $cos_winner_lng;
$north_y = -$local_grid_step * $sin_winner_lat * $sin_winner_lng;
$north_z = $local_grid_step * $cos_winner_lat;
if ($local_grid_radius > 0 && $local_grid_step > 0) {
$r = intval($local_grid_radius / $local_grid_step);
$r2 = $r * $r;
$h = $local_kernel_radius;
$h2 = $h * $h;
$max_L_sum = -1;
$max_L_sum_idx = -1;
for ($i = -$r; $i <= $r; $i++) {
$winner_east_x = $winner_x + $i * $east_x;
$winner_east_y = $winner_y + $i * $east_y;
$winner_east_z = $winner_z + $i * $east_z;
$j_max = intval(sqrt($r2 - $i * $i));
for ($j = -$j_max; $j <= $j_max; $j++) {
$x = $winner_east_x + $j * $north_x;
$y = $winner_east_y + $j * $north_y;
$z = $winner_east_z + $j * $north_z;
$L_sum = 0;
for ($k = 0; $k < $lat_lng_count; $k++) {
$dx = $x - $x_arr[$k];
$dy = $y - $y_arr[$k];
$dz = $z - $z_arr[$k];
$d2 = $dx * $dx + $dy * $dy + $dz * $dz;
if ($d2 < $h2) {
$L_sum += $w_arr[$k] * ($h2 - $d2); // Epanechnikov
}
}
if ($max_L_sum < $L_sum) {
$max_L_sum = $L_sum;
$max_L_sum_i = $i;
$max_L_sum_j = $j;
}
}
}
$x = $winner_x + $max_L_sum_i * $east_x + $max_L_sum_j * $north_x;
$y = $winner_y + $max_L_sum_i * $east_y + $max_L_sum_j * $north_y;
$z = $winner_z + $max_L_sum_i * $east_z + $max_L_sum_j * $north_z;
} else if ($local_grid_radius > 0) {
$r = $local_grid_radius;
$r2 = $r * $r;
$wx_sum = 0;
$wy_sum = 0;
$wz_sum = 0;
$w_sum = 0;
for ($k = 0; $k < $lat_lng_count; $k++) {
$xk = $x_arr[$k];
$yk = $y_arr[$k];
$zk = $z_arr[$k];
$dx = $winner_x - $xk;
$dy = $winner_y - $yk;
$dz = $winner_z - $zk;
$d2 = $dx * $dx + $dy * $dy + $dz * $dz;
if ($d2 <= $r2) {
$wk = $w_arr[$k];
$wx_sum += $wk * $xk;
$wy_sum += $wk * $yk;
$wz_sum += $wk * $zk;
$w_sum += $wk;
}
}
$x = $wx_sum / $w_sum;
$y = $wy_sum / $w_sum;
$z = $wz_sum / $w_sum;
$max_L_sum_i = false;
$max_L_sum_j = false;
} else {
return array($winner_lat, $winner_lng, $max_K_sum_i, false, false);
}
$deg = 180 / M_PI;
$a2 = $a * $a;
$e4 = $e2 * $e2;
$p = sqrt($x * $x + $y * $y);
$zeta = (1 - $e2) * $z * $z / $a2;
$rho = ($p * $p / $a2 + $zeta - $e4) / 6;
$rho3 = $rho * $rho * $rho;
$s = $e4 * $zeta * $p * $p / (4 * $a2);
$t = pow($s + $rho3 + sqrt($s * ($s + 2 * $rho3)), 1 / 3);
$u = $rho + $t + $rho * $rho / $t;
$v = sqrt($u * $u + $e4 * $zeta);
$w = $e2 * ($u + $v - $zeta) / (2 * $v);
$k = 1 + $e2 * (sqrt($u + $v + $w * $w) + $w) / ($u + $v);
$lat = atan($k * $z / $p) * $deg;
$lng = atan2($y, $x) * $deg;
return array($lat, $lng, $max_K_sum_i, $max_L_sum_i, $max_L_sum_j);
}
距离是欧几里德而不是大圆的事实应该对手头的任务产生微不足道的影响.计算大圆距离会更加麻烦,并且只会导致非常远的点的重量显着降低 - 但这些点的重量已经非常低.原则上,不同的内核可以实现相同的效果.具有超出一定距离的完全截止的内核,如Epanechnikov内核,根本没有这个问题(在实践中).
WGS84基准的lat,lng和x,y,z之间的转换(尽管不保证数值稳定性)更多地作为参考,而不是真正的需要.如果要考虑高度,或者需要更快的反向转换,请参阅维基百科文章.
Epanechnikov内核除了比高斯和Trossian内核"更本地"外,还具有第二个循环最快的优点,即O(ng),其中g是局部网格的点数,并且可以如果n很大,也可以在第一个循环中使用,即O(n²).
找到密度估计值最大的点.
应该非常简单明了.使用大致覆盖直径大型机场的核心半径.2D Gaussian或Epanechnikov内核应该没问题.
http://en.wikipedia.org/wiki/Multivariate_kernel_density_estimation
这类似于计算堆映射:http://en.wikipedia.org/wiki/Heat_map ,然后在那里找到最亮点.除了它立即计算亮度.
为了好玩,我将DBpedia的Geocoordinates(即维基百科)的1%样本读入ELKI,将其投影到3D空间并启用密度估计叠加(隐藏在可视化器散点图菜单中).你可以看到欧洲有一个热点,而在美国则有一个较小的热点.我相信欧洲的热点是波兰.最后我查了一下,显然有人在波兰的几乎任何一个城镇为Geocoordinates创建了一篇维基百科文章.遗憾的是,ELKI可视化工具既不允许您放大,旋转或减少内核带宽,也无法直观地找到最密集的点.但实施自己很简单; 你可能也不需要进入3D空间,但可以只使用纬度和经度.
核密度估计应该在大量应用中可用.R中的那个可能更强大.我刚刚在ELKI中发现了这个热图,所以我知道如何快速访问它.有关相关的R函数,请参见http://stat.ethz.ch/R-manual/R-devel/library/stats/html/density.html.
在您的数据中,在R中,尝试例如:
library(kernSmooth)
smoothScatter(data, nbin=512, bandwidth=c(.25,.25))
这应该表明对芝加哥的强烈偏好.
library(kernSmooth)
dens=bkde2D(data, gridsize=c(512, 512), bandwidth=c(.25,.25))
contour(dens$x1, dens$x2, dens$fhat)
maxpos = which(dens$fhat == max(dens$fhat), arr.ind=TRUE)
c(dens$x1[maxpos[1]], dens$x2[maxpos[2]])
收益率[1] 42.14697 -88.09508,距离芝加哥机场不到10英里.
为了获得更好的坐标尝试:
- 在估计坐标周围的20x20英里区域重新运行
- 该区域内的非分箱KDE
- 更好的带宽选择
dpik - 更高的网格分辨率
在天体物理学中,我们使用所谓的"半质量半径".给定分布及其中心,半质量半径是包含分布的一半点的圆的最小半径.
该数量是点分布的特征长度.
如果你想要直升机的家是最大点集中的地方,那么它是具有最小半质量半径的点!
我的算法如下:对于每个点,您计算这个半质量半径,以当前点的分布为中心.直升机的"家"将是具有最小半质量半径的点.
我已经实现了它并且计算中心42.149994 -88.133698(在芝加哥)我也使用了总质量的0.2而不是通常用于天体物理学的0.5(一半).
这是我(在python中)alghorithm找到直升飞机的家:
import math
import numpy
def inside(points,center,radius):
ids=(((points[:,0]-center[0])**2.+(points[:,1]-center[1])**2.)<=radius**2.)
return points[ids]
points = numpy.loadtxt(open('points.txt'),comments='#')
npoints=len(points)
deltar=0.1
idcenter=None
halfrmin=None
for i in xrange(0,npoints):
center=points[i]
radius=0.
stayHere=True
while stayHere:
radius=radius+deltar
ninside=len(inside(points,center,radius))
#print 'point',i,'r',radius,'in',ninside,'center',center
if(ninside>=npoints*0.2):
if(halfrmin==None or radius<halfrmin):
halfrmin=radius
idcenter=i
print 'point',i,halfrmin,idcenter,points[idcenter]
stayHere=False
#print halfrmin,idcenter
print points[idcenter]