Как найти центр кластера точек данных?

Ответ 1

Следующее решение работает, даже если точки рассеяны по всей Земле, путем преобразования широты и долготы в декартовы координаты. Он делает своего рода KDE (оценка плотности ядра), но в первом проходе сумма ядер оценивается только в точках данных. Ядро должно быть выбрано для соответствия этой проблеме. В приведенном ниже коде это то, что я мог в шутку/самонадеянно назвать тросианским, т.е. 2-d²/h² для d≤h и h²/d² для d > h (где d - евклидово расстояние, h - "ширина полосы" $global_kernel_radius), но он также может быть гауссовым (e ^-d²/2h²), ядром Эпанечникова (1-d²/h² для d < h, 0 в противном случае) или другим ядром. Дополнительный второй проход уточняет поиск локально, либо путем суммирования независимого ядра на локальной сетке, либо путем вычисления центроида в обоих случаях в окружении, определяемом $local_grid_radius.

В сущности каждая точка суммирует все точки вокруг (включая себя), взвешивая их больше, если они ближе (по кривой колокола), а также взвешивает их с помощью дополнительной весовой матрицы $w_arr. Победителем является точка с максимальной суммой. Как только победитель будет найден, "дом", который мы ищем, можно найти, повторяя тот же процесс локально вокруг победителя (используя другую колоколообразную кривую), или его можно оценить как "центр масс" всех точек в пределах данного радиуса от победителя, где радиус может быть равен нулю.

Алгоритм должен быть адаптирован к задаче, выбирая подходящие ядра, выбирая, как уточнить поиск локально, и путем настройки параметров. Для примера набора данных может быть хорошей отправной точкой ядро ​​Троссана для первого прохода и ядро ​​Эпанечникова для второго прохода со всеми тремя радиусами, установленными в 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);
}

Тот факт, что расстояния являются евклидовыми, а не большим кругом, должен иметь незначительные эффекты для этой задачи. Вычисление расстояний с большим кругом было бы намного более громоздким и приводило бы к значительному снижению веса очень далеких точек, но эти очки уже имеют очень низкий вес. В принципе, тот же эффект может быть достигнут другим ядром. Ядра, которые имеют полное отсечение на некотором расстоянии, например ядро ​​Эпанечникова, вообще не имеют этой проблемы (на практике).

Преобразование между lat, lng и x, y, z для базы данных WGS84 дается точно (хотя и без гарантии численной устойчивости) скорее как ссылка, чем из-за истинной потребности. Если высота должна быть принята во внимание или требуется более быстрое обратное преобразование, обратитесь к статье Wikipedia.

Ядро Эпанечникова, помимо "более локального", чем ядра Гаусса и Троссиана, имеет то преимущество, что оно является самым быстрым для второго цикла, который является O (ng), где g - число точек локальной сетки, и может также использоваться в первом цикле, который является O (n²), если n велико.

Ответ 2

Это можно решить, найдя под угрозу поверхность. См. Формула Rossmo.

Это проблема хищника. Учитывая набор географически расположенных туш, где находится логово хищника? Формула Rossmo решает эту проблему.

Ответ 3

Найдите точку с оценкой наибольшей плотности.

Должно быть довольно просто. Используйте радиус ядра, который примерно покрывает большой аэропорт в диаметре. Ядро 2D Gaussian или Epanechnikov должно быть в порядке.

http://en.wikipedia.org/wiki/Multivariate_kernel_density_estimation

Это похоже на вычисление карты кучи: http://en.wikipedia.org/wiki/Heat_map а затем найти там самое яркое место. За исключением того, что он сразу же вычисляет яркость.

Для удовольствия я прочитал 1% -ный образец геокоординатов DBpedia (т.е. Wikipedia) в ELKI, проецировал его в 3D-пространство и включил наложение оценки плотности (скрытое в меню рассеивателя визуализаторов). Вы видите, что в Европе есть горячая точка, и в меньшей степени в США. Я считаю, что горячая точка в Европе - это Польша. Последнее, что я проверил, кто-то, видимо, создал статью в Википедии с Geocoordinates для почти любого города в Польше. К сожалению, визуализатор ELKI не позволяет увеличить, повернуть или уменьшить пропускную способность ядра, чтобы визуально найти самую плотную точку. Но это просто осуществить сами; вам, вероятно, также не нужно идти в 3D-пространство, но может просто использовать широты и долготы.

Оценка плотности ядра должна быть доступна в тоннах приложений. Вероятно, один из R является гораздо более мощным. Недавно я обнаружил эту тепловую карту в ELKI, поэтому я знал, как быстро ее получить. См. http://stat.ethz.ch/R-manual/R-devel/library/stats/html/density.html для связанной функции R.

В ваших данных, в 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
• более высокое разрешение сетки

Ответ 4

в астрофизике мы используем так называемый радиус полумассы. Учитывая распределение и его центр, радиус полумассы - это минимальный радиус круга, который содержит половину точек вашего распределения.

Эта величина является характерной длиной распределения точек.

Если вы хотите, чтобы дом вертолета находился там, где точки максимально сконцентрированы, так что это точка, которая имеет минимальный радиус полумассы!

Мой алгоритм выглядит следующим образом: для каждой точки вы вычисляете этот радиус половинной массы, центрирующий распределение в текущей точке. "Дом" вертолета будет точкой с минимальным радиусом половинной массы.

Я реализовал его, и вычисленный центр 42.149994 -88.133698 (который находится в Чикаго) Я также использовал 0,2 общей массы вместо 0,5 (половины), обычно используемых в астрофизике.

Это мой (в python) альгорифм, который находит дом вертолета:

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]

Ответ 5

Вы можете использовать DBSCAN для этой задачи.

DBSCAN - это кластеризация на основе плотности с понятием шума. Вам нужны два параметра:

Сначала количество точек, которые кластер должен иметь минимум "minpoints". И второй параметр окрестности, называемый "epsilon", который устанавливает порог расстояния к окружающим точкам, которые должны быть включены в ваш кластер.

Весь алгоритм работает следующим образом:

• Начните с произвольной точки в вашем наборе, которая еще не была посещена.
• Получить точки из окрестности эпсилона отметьте все, как посетили
  • если вы нашли достаточно точек в этой окрестности ( > параметр minpoints), вы запускаете новый кластер и назначаете эти точки. Теперь снова переходите к шагу 2 для каждой точки в этом кластере.
  • если у вас его нет, объявите этот пункт как шум
• перейдите снова, пока вы не посетили все точки.

Это действительно просто реализовать, и есть много фреймворков, которые уже поддерживают этот алгоритм. Чтобы найти среднее значение вашего кластера, вы можете просто взять среднее значение всех назначенных точек из его окрестности.

Однако, в отличие от метода, который предлагает @TylerDurden, для этого требуется параметризация, поэтому вам нужно найти некоторые настроенные вручную параметры, которые соответствуют вашей проблеме.

В вашем случае вы можете попытаться установить минимальные точки на 10% от ваших общих очков, если самолет, скорее всего, останется в 10% от времени, когда вы проедете в аэропорту. Параметр плотности epsilon зависит от разрешения вашего географического датчика и используемой метрики расстояния. Я бы предложил расстояние haversine для географических данных.

Ответ 6

Все, что у меня на этой машине, это старый компилятор, поэтому я сделал версию ASCII. Он "рисует" (в ASCII) точки карты - это точки, X - это источник реального источника, где G - это предполагаемый источник. Если два перекрытия, отображается только X.

Примеры (ТРУДНОСТЬ 1.5 и 3 соответственно):

Точки генерируются путем выбора случайной точки в качестве источника, затем случайного распределения точек, что делает их более вероятными ближе к источнику.

DIFFICULTY - постоянная с плавающей запятой, которая регулирует генерацию исходной точки - насколько вероятнее, что точки должны быть ближе к источнику - если это 1 или меньше, программа должна уметь точно определять источник, или очень близко. В 2.5 он все равно будет довольно приличным. В 4+ он начнет угадываться хуже, но я думаю, что он все еще догадывается лучше, чем человек.

Его можно оптимизировать, используя бинарный поиск по X, тогда Y - это сделало бы предположение хуже, но было бы намного, намного быстрее. Или, начиная с более крупных блоков, затем разбивая лучший блок дальше (или лучший блок и 8, окружающие его). Для системы с более высоким разрешением потребуется одна из них. Однако это довольно наивный подход, но, похоже, он хорошо работает в системе 80x24.: D

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define Y 24
#define X 80

#define DIFFICULTY 1 // Try different values... 

static int point[Y][X];

double dist(int x1, int y1, int x2, int y2)
{
    return sqrt((y1 - y2)*(y1 - y2) + (x1 - x2)*(x1 - x2));
}

main()
{
    srand(time(0));
    int y = rand()%Y;
    int x = rand()%X;

    // Generate points
    for (int i = 0; i < Y; i++)
    {
        for (int j = 0; j < X; j++)
        {
            double u = DIFFICULTY * pow(dist(x, y, j, i), 1.0 / DIFFICULTY);
            if ((int)u == 0)
                u = 1;
            point[i][j] = !(rand()%(int)u);
        }
    }

    // Find best source
    int maxX = -1;
    int maxY = -1;
    double maxScore = -1;
    for (int cy = 0; cy < Y; cy++)
    {
        for (int cx = 0; cx < X; cx++)
        {
            double score = 0;
            for (int i = 0; i < Y; i++)
            {
                for (int j = 0; j < X; j++)
                {
                    if (point[i][j] == 1)
                    {
                        double d = dist(cx, cy, j, i);
                        if (d == 0)
                            d = 0.5;
                        score += 1000 / d;
                    }
                }
            }
            if (score > maxScore || maxScore == -1)
            {
                maxScore = score;
                maxX = cx;
                maxY = cy;
            }
        }
    }

    // Print out results
    for (int i = 0; i < Y; i++)
    {
        for (int j = 0; j < X; j++)
        {
            if (i == y && j == x)
                printf("X");
            else if (i == maxY && j == maxX)
                printf("G");            
            else if (point[i][j] == 0)
                printf(" ");
            else if (point[i][j] == 1)
                printf(".");
        }
    }
    printf("Distance from real source: %f", dist(maxX, maxY, x, y));

    scanf("%d", 0);

} 

Ответ 7

Как разделить карту на многие зоны, а затем найти центр плоскости в зоне с самой плоской. Алгоритм будет что-то вроде этого

set Zones[40]
foreach Plane in Planes
   Zones[GetZone(Plane.position)].Add(Plane)

set MaxZone = Zones[0]
foreach Zone in Zones
   if MaxZone.Length() < Zone.Length()
           MaxZone = Zone

set Center
foreach Plane in MaxZone
     Center.X += Plane.X
     Center.Y += Plane.Y
Center.X /= MaxZone.Length
Center.Y /= MaxZone.Length

Ответ 8

Виртуальная земля имеет очень хорошее объяснение того, как вы можете сделать это относительно быстро. Они также предоставили примеры кода. Пожалуйста, посмотрите http://soulsolutions.com.au/Articles/ClusteringVirtualEarthPart1.aspx

Ответ 9

Простая модель смеси, похоже, очень хорошо справляется с этой проблемой.

В общем, чтобы получить точку, которая минимизирует расстояние до всех других точек в наборе данных, вы можете просто взять среднее значение. В этом случае вы хотите найти точку, которая минимизирует расстояние от подмножества сосредоточенных точек. Если вы постулируете, что точка может быть либо из концентрированного набора точек интереса, либо из диффузного набора фоновых точек, то это дает модель смеси.

Я включил код python ниже. Концентрированная область моделируется высокоточным нормальным распределением, а фоновая точка моделируется либо низкочастотным нормальным распределением, либо равномерным распределением по ограничивающей рамке на наборе данных (имеется строка кода, которая может быть прокомментирована переключение между этими параметрами). Кроме того, модели смеси могут быть несколько неустойчивыми, поэтому несколько раз запускайте алгоритм ЭМ с случайными начальными условиями и выбирая прогон с наивысшим лог-правдоподобием, дают лучшие результаты.

Если вы на самом деле смотрите на самолеты, то добавление какой-то зависящей от времени динамики, вероятно, улучшит вашу способность очень глубоко обосновать базовую базу.

Я также опасаюсь формулы Россимо, потому что в нее входят некоторые довольно сильные предположения о распределении преступлений.

#the dataset
sdata='''41.892694,-87.670898
42.056048,-88.000488
41.941744,-88.000488
42.072361,-88.209229
42.091933,-87.982635
42.149994,-88.133698
42.171371,-88.286133
42.23241,-88.305359
42.196811,-88.099365
42.189689,-88.188629
42.17646,-88.173523
42.180531,-88.209229
42.18168,-88.187943
42.185496,-88.166656
42.170485,-88.150864
42.150634,-88.140564
42.156743,-88.123741
42.118555,-88.105545
42.121356,-88.112755
42.115499,-88.102112
42.119319,-88.112411
42.118046,-88.110695
42.117791,-88.109322
42.182189,-88.182449
42.194145,-88.183823
42.189057,-88.196182
42.186513,-88.200645
42.180917,-88.197899
42.178881,-88.192062
41.881656,-87.6297
41.875521,-87.6297
41.87872,-87.636566
41.872073,-87.62661
41.868239,-87.634506
41.86875,-87.624893
41.883065,-87.62352
41.881021,-87.619743
41.879998,-87.620087
41.8915,-87.633476
41.875163,-87.620773
41.879125,-87.62558
41.862763,-87.608757
41.858672,-87.607899
41.865192,-87.615795
41.87005,-87.62043
42.073061,-87.973022
42.317241,-88.187256
42.272546,-88.088379
42.244086,-87.890625
42.044512,-88.28064
39.754977,-86.154785
39.754977,-89.648437
41.043369,-85.12207
43.050074,-89.406738
43.082179,-87.912598
42.7281,-84.572754
39.974226,-83.056641
38.888093,-77.01416
39.923692,-75.168457
40.794318,-73.959961
40.877439,-73.146973
40.611086,-73.740234
40.627764,-73.234863
41.784881,-71.367187
42.371988,-70.993652
35.224587,-80.793457
36.753465,-76.069336
39.263361,-76.530762
25.737127,-80.222168
26.644083,-81.958008
30.50223,-87.275391
29.436309,-98.525391
30.217839,-97.844238
29.742023,-95.361328
31.500409,-97.163086
32.691688,-96.877441
32.691688,-97.404785
35.095754,-106.655273
33.425138,-112.104492
32.873244,-117.114258
33.973545,-118.256836
33.681497,-117.905273
33.622982,-117.734985
33.741828,-118.092041
33.64585,-117.861328
33.700707,-118.015137
33.801189,-118.251343
33.513132,-117.740479
32.777244,-117.235107
32.707939,-117.158203
32.703317,-117.268066
32.610821,-117.075806
34.419726,-119.701538
37.750358,-122.431641
37.50673,-122.387695
37.174817,-121.904297
37.157307,-122.321777
37.271492,-122.033386
37.435238,-122.217407
37.687794,-122.415161
37.542025,-122.299805
37.609506,-122.398682
37.544203,-122.0224
37.422151,-122.13501
37.395971,-122.080078
45.485651,-122.739258
47.719463,-122.255859
47.303913,-122.607422
45.176713,-122.167969
39.566,-104.985352
39.124201,-94.614258
35.454518,-97.426758
38.473482,-90.175781
45.021612,-93.251953
42.417881,-83.056641
41.371141,-81.782227
33.791132,-84.331055
30.252543,-90.439453
37.421248,-122.174835
37.47794,-122.181702
37.510628,-122.254486
37.56943,-122.346497
37.593373,-122.384949
37.620571,-122.489319
36.984249,-122.03064
36.553017,-121.893311
36.654442,-121.772461
36.482381,-121.876831
36.15042,-121.651611
36.274518,-121.838379
37.817717,-119.569702
39.31657,-120.140991
38.933041,-119.992676
39.13785,-119.778442
39.108019,-120.239868
38.586082,-121.503296
38.723354,-121.289062
37.878444,-119.437866
37.782994,-119.470825
37.973771,-119.685059
39.001377,-120.17395
40.709076,-73.948975
40.846346,-73.861084
40.780452,-73.959961
40.778829,-73.958931
40.78372,-73.966012
40.783688,-73.965325
40.783692,-73.965615
40.783675,-73.965741
40.783835,-73.965873
'''

import StringIO
import numpy as np
import re

import matplotlib.pyplot as plt

def lp(l):
    return map(lambda m: float(m.group()),re.finditer('[^, \n]+',l))

data=np.array(map(lp,StringIO.StringIO(sdata)))

xmn=np.min(data[:,0])
xmx=np.max(data[:,0])
ymn=np.min(data[:,1])
ymx=np.max(data[:,1])

# area of the point set bounding box
area=(xmx-xmn)*(ymx-ymn)

M_ITER=100 #maximum number of iterations
THRESH=1e-10 # stopping threshold

def em(x):
    print '\nSTART EM'
    mlst=[]

    mu0=np.mean( data , 0 ) # the sample mean of the data - use this as the mean of the low-precision gaussian

    # the mean of the high-precision Gaussian - this is what we are looking for
    mu=np.random.rand( 2 )*np.array([xmx-xmn,ymx-ymn])+np.array([xmn,ymn])

    lam_lo=.001  # precision of the low-precision Gaussian
    lam_hi=.1 # precision of the high-precision Gaussian
    prz=np.random.rand( 1 ) # probability of choosing the high-precision Gaussian mixture component

    for i in xrange(M_ITER):
        mlst.append(mu[:])

        l_hi=np.log(prz)+np.log(lam_hi)-.5*lam_hi*np.sum((x-mu)**2,1)
        #low-precision normal background distribution
        l_lo=np.log(1.0-prz)+np.log(lam_lo)-.5*lam_lo*np.sum((x-mu0)**2,1)
        #uncomment for the uniform background distribution
        #l_lo=np.log(1.0-prz)-np.log(area)

        #expectation step
        zs=1.0/(1.0+np.exp(l_lo-l_hi))

        #compute bound on the likelihood 
        lh=np.sum(zs*l_hi+(1.0-zs)*l_lo)
        print i,lh

        #maximization step
        mu=np.sum(zs[:,None]*x,0)/np.sum(zs) #mean
        lam_hi=np.sum(zs)/np.sum(zs*.5*np.sum((x-mu)**2,1)) #precision
        prz=1.0/(1.0+np.sum(1.0-zs)/np.sum(zs)) #mixure component probability

        try:
            if np.abs((lh-old_lh)/lh)<THRESH:
                break
        except: 
            pass

        old_lh=lh

        mlst.append(mu[:])

    return lh,lam_hi,mlst    

if __name__=='__main__':

    #repeat the EM algorithm a number of times and get the run with the best log likelihood
    mx_prm=em(data)
    for i in xrange(4):
        prm=em(data)

        if prm[0]>mx_prm[0]:
            mx_prm=prm

        print prm[0]
        print mx_prm[0]

    lh,lam_hi,mlst=mx_prm
    mu=mlst[-1]

    print 'best loglikelihood:', lh
    #print 'final precision value:', lam_hi
    print 'point of interest:', mu
    plt.plot(data[:,0],data[:,1],'.b')

    for m in mlst:
        plt.plot(m[0],m[1],'xr')

    plt.show()

Ответ 10

Вы можете легко приспособить формулу Rossmo, процитированную Tyler Durden, к вашему делу с несколькими простыми примечаниями:

Формула:

Эта формула дает что-то близкое к вероятности наличия базовой операции для хищника или серийного убийцы. В вашем случае это может дать вероятность того, что база будет в определенной точке. Я объясню позже, как его использовать. U может записать это так:

Proba (основание в точке A) = Sum {во всех точках} (Phi/(dist ^ f) + (1-Phi) (B * (gf))/(2B-dist) ^ g)

Использование евклидова расстояния

Вы хотите, чтобы евклидское расстояние, а не манхэттенское, потому что самолёт или вертолет не связаны с дорогами/улицами. Таким образом, использование евклидовой дистанции является правильным способом, если вы отслеживаете самолет, а не серийный убийца. Таким образом, "dist" в формуле является эвклидовым расстоянием между точечным ур-тестированием и рассматриваемым пятном

Принимая разумную переменную B

Переменная B использовалась для представления правила: "разумно умный убийца не убивает своего соседа". В вашем случае будет применяться также, потому что никто не использует самолет /roflcopter, чтобы добраться до следующего угла улицы. мы можем предположить, что минимальное путешествие, например, 10 км или что-нибудь разумное при применении к вашему делу.

Экспоненциальный коэффициент f

Фактор f используется для добавления веса к расстоянию. Например, если все пятна находятся в небольшой области, вам может понадобиться большой коэффициент f, так как вероятность того, что аэропорт/база/HQ будет быстро уменьшаться, если все ваши данные будут находиться в одном секторе. g работает аналогичным образом, позволяет выбрать размер "база вряд ли будет рядом с областью"

Фактор Phi:

Снова этот фактор необходимо определить, используя свои знания о проблеме. Это позволяет выбрать наиболее точный коэффициент между "базой близок к пятнам" и "я не буду использовать плоскость, чтобы сделать 5 м", если, например, вы думаете, что вторая почти не имеет значения, вы можете установить Phi на 0,95 (0<Phi<1) Если оба интересны, phi будет около 0,5

Как реализовать его как нечто полезное:

Сначала вы хотите разделить свою карту на маленькие квадраты: завязать карту (точно так же, как и invisal) (чем меньше квадраты, тем точнее результат (в общем)), а затем используя формулу, чтобы найти более вероятное местоположение. На самом деле сетка - это всего лишь массив со всеми возможными местоположениями. (если вы хотите быть точным, вы увеличиваете количество возможных пятен, но для этого потребуется больше вычислительного времени, а PhP не известен своей удивительной скоростью).

Алгоритм:

//define all the factors you need(B , f , g , phi)

for(i=0..mesh_size) // computing the probability of presence for each square of the mesh
{
  P(i)=0;
  geocode squarePosition;//GeoCode of the square center 
  for(j=0..geocodearray_size)//sum on all the known spots
  {
     dist=Distance(geocodearray[j],squarePosition);//small function returning distance between two geocodes

         P(i)+=(Phi/pow(dist,f))+(1-Phi)*pow(B,g-f)/pow(2B-dist,g);
  }
 }

return geocode corresponding to max(P(i))

Надеюсь, что это поможет вам

Ответ 11

Сначала я хотел бы выразить свою любовь к вашему методу для иллюстрации и объяснения проблемы.

Если бы я был на твоем месте, я бы выбрал алгоритм плотности плотности, например DBSCAN а затем после кластеризации областей и удаления точек шума останется несколько областей (вариантов). Тогда я возьму кластер с наивысшей плотностью точек и вычислить среднюю точку и найти ближайшую действительную точку. сделал, нашел это место!:.)

Привет,

Ответ 12

Почему бы не что-то вроде этого:

• Для каждой точки вычислите расстояние до всех остальных точек и суммируйте их.
• Точка с наименьшей суммой - это ваш центр.

Может быть, сумма не самая лучшая метрика для использования. Возможно, точка с самыми "небольшими расстояниями"?

Ответ 13

Сумма на расстояния. Возьмите точку с наименьшим суммарным расстоянием.

function () {
    for i in points P:
        S[i] = 0
        for j in points P:
            S[i] += distance(P[i], P[j])
    return min(S);
}

Ответ 14

Вы можете взять минимальное связующее дерево и удалить самые длинные края. Меньшие деревья дают вам центр для поиска. Имя алгоритма представляет собой односвязную k-кластеризацию. Здесь есть сообщение: https://stats.stackexchange.com/info/1475/visualization-software-for-clustering.