Wie finde ich die Mitte aus einem cluster von Daten-Punkten?

Folgende Lösung funktioniert auch, wenn die Punkte sind überall auf der Erde, durch die Umwandlung von Längen-und Breitengrad zu kartesischen Koordinaten. Es funktioniert eine Art von KDE (kernel density estimation), aber in einem ersten Durchlauf die Summe der Kerne ausgewertet wird nur die Daten-Punkte. Der kernel sollte so gewählt werden, passen Sie das problem. Im code unten ist es das, was ich mir scherzhaft/überheblich nennen eine Trossian, d.h., 2-d2/h2 für d≤h und h2/d2 für d>h (wobei d der euklidische Abstand und h ist die "Bandbreite" $global_kernel_radius), aber es könnte auch eine Gauß - (e^-d2/2h2), ein Epanechnikov-kernel (1-d2/h2 für d<h, 0 sonst), oder einen anderen kernel. Eine optionale zweite Durchlauf verfeinert die Suche, lokal, entweder durch addieren einer unabhängigen kernel auf einem lokalen Netz, oder durch die Berechnung der Schwerpunkt, in beiden Fällen in eine Umgebung definiert, die von $local_grid_radius.

Im wesentlichen, jeder Punkt Summen all die Punkte, die es hat, um mit sich selbst, wiegt Sie mehr, wenn Sie näher (durch die bell-Kurve), und auch Wiegen, die Ihnen durch die optional Gewicht array $w_arr. Der Gewinner ist der Punkt mit der maximalen Summe. Nachdem der Sieger gefunden wurde, wird die "home" suchen wir für die gefunden werden können, wiederholen Sie den gleichen Prozess lokal um den Gewinner (mit einem anderen bell-Kurve), oder es kann geschätzt werden, werden die "center of mass" aller Punkte innerhalb eines bestimmten radius um den Gewinner, wobei der radius null werden kann.

Muss der Algorithmus angepasst werden, um das problem, indem Sie die entsprechenden Kernel, indem Sie, wie um die Suche zu verfeinern lokal, und durch einstellen der Parameter. Für das Beispiel-dataset, das Trossian kernel für den ersten pass, und der Epanechnikov-kernel für den zweiten pass, alle 3 Radien auf 30 mi und grid step 1 mi könnte ein guter Ausgangspunkt, aber nur, wenn die beiden sub-Clustern von Chicago sollte gesehen werden als eine große Gruppe. Ansonsten kleinere Radien müssen so gewählt werden.

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);
}

Die Tatsache, dass Entfernungen entsprechen der euklidischen und nicht groß-Kreis sollte nur vernachlässigbare Auswirkungen für die Aufgabe zur hand. Die Berechnung der great-circle Entfernungen wäre viel umständlicher und würde nur dazu führen das Gewicht sehr weit Punkte werden deutlich geringer - aber diese Punkte haben schon einen sehr niedrigen Gewicht. Im Prinzip der gleiche Effekt könnte erzielt werden, indem ein anderer kernel. Kerne, die haben eine komplette cut-off-jenseits einiger Entfernung, wie der Epanechnikov-kernel, haben dieses problem nicht bei allen (in der Praxis).

Die Umrechnung zwischen lat,lng und x,y,z für das WGS84 datum exakt angegeben (allerdings ohne Gewährleistung der numerischen Stabilität) mehr als eine Referenz ist, als durch eine echte Notwendigkeit. Wenn die Höhe berücksichtigt werden, oder wenn eine schnellere back-Konvertierung erforderlich ist, entnehmen Sie bitte der Wikipedia-Artikel.

Den Epanechnikov-kernel, abgesehen davon, dass "mehr lokales" als die Gauß-und Trossian-Kernel, hat den Vorteil, dass das Schnellste für die zweite Schleife ist O(ng), wobei g die Anzahl der Punkte, die von dem lokalen Netz, und kann auch eingesetzt werden, in die erste Schleife ist O(n2), wenn n groß ist.

Finden Sie den Punkt mit der größte Dichte Schätzung.

Sein sollte ziemlich einfach. Verwenden Sie einen kernel-radius, der etwa deckt einen großen Flughafen im Durchmesser. Einen 2D-GAUSS-oder Epanechnikov-kernel sollte in Ordnung sein.

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

Dies ist ähnlich wie die Berechnung einer Heap-Karte: http://en.wikipedia.org/wiki/Heat_map
und dann das finden der hellsten Stelle. Außer es berechnet die Helligkeit sofort.

For fun habe ich gelesen, eine 1% - Stichprobe der geografischen Koordinaten von DBpedia (z.B. Wikipedia) in ELKI, projiziert es in den 3D-Raum und aktiviert die Dichte Schätzung overlay (versteckt in den Visualisierungen scatterplot-Menü). Sie können sehen, es ist ein hotspot in Europa, und in geringerem Ausmaß in den USA. Der hotspot in Europa ist Polen, glaube ich. Zuletzt habe ich geprüft, jemand hatte offenbar erstellt einen Wikipedia-Artikel mit geografischen Koordinaten für so ziemlich jede Stadt in Polen. Die ELKI-visualizer, leider, weder können Sie Zoomen, drehen, oder reduzieren Sie die kernel-Bandbreite optisch finden die meisten dichten Punkt. Aber es ist einfach zu implementieren Sie sich selbst; Sie wahrscheinlich auch nicht brauchen, um in den 3D-Raum, kann aber nur verwenden breiten-und Längengrade.

Kernel-Dichte-Schätzung sollte in Tonnen von Anwendungen. Die eins in R ist wahrscheinlich viel stärker. Ich habe vor kurzem entdeckt, diese heatmap im ELKI, so dass ich wusste, wie schnell darauf zugreifen. Siehe z.B. http://stat.ethz.ch/R-manual/R-devel/library/stats/html/density.html für einen zugehörigen R-Funktion.

Auf Ihre Daten, R, versuchen Sie es zum Beispiel:

library(kernSmooth)
smoothScatter(data, nbin=512, bandwidth=c(.25,.25))

dies sollte zeigen eine starke Präferenz für Chicago.

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]])

Erträge [1] 42.14697 -88.09508 - das ist weniger als 10 Meilen vom Flughafen von Chicago.

Besser koordiniert versuchen:

• rerunning auf einer 20x20 km-Bereich um den geschätzten Koordinaten
• eine non-binned KDE in diesem Bereich
• bessere Bandbreite mit dpik
• höhere Auflösung

Können Sie DBSCAN für diese Aufgabe.

DBSCAN ist ein Dichte-basiertes clustering mit einer Vorstellung von Lärm. Sie benötigen zwei Parameter:

Zunächst die Anzahl der Punkte zu einem cluster sollte mindestens "minpoints".
Und zweitens, ein Nachbarschafts-parameter genannt "epsilon" setzt einen Abstand Schwelle zu den umliegenden Punkte, die enthalten sein sollten in Ihrem cluster.

Den gesamten Algorithmus funktioniert wie folgt:

1. Starten mit einem beliebigen Punkt in Ihrem set, der bisher nicht besucht haben, noch
2. Abrufen Punkte aus dem epsilon Nachbarschaft markieren Sie alle besucht
   1. wenn Sie gefunden haben, genug Punkte in diesem Viertel (> minpoints parameter) starten Sie einen neuen cluster und weisen Sie diese Punkte. Jetzt recurse in Schritt 2 wieder für jeden Punkt in diesem cluster.
   2. wenn Sie noch nicht besitzen, erklären diesen Punkt als Rauschen
3. gehen alle immer wieder, bis Sie besucht haben, alle Punkte

Es ist wirklich einfach umzusetzen und es gibt viele frameworks, die Unterstützung dieser Algorithmus bereits. Zu finden, den Mittelwert des Clusters sind, können Sie nehmen Sie einfach der Mittelwert aller zugeordneten Punkte aus Ihrer Nachbarschaft.

Jedoch, im Gegensatz zu der Methode, die @TylerDurden schlägt vor, dies braucht eine Parametrisierung - also müssen Sie etwas finden, von hand gestimmt Parameter passend zu Ihrem problem.

In Ihrem Fall können Sie versuchen die minpoints zu 10% eurer gesamten Punkte, wenn das Flugzeug wird voraussichtlich bleiben 10% der Zeit, die Sie verfolgen auf einem Flughafen. Die Dichte-parameter epsilon, hängt von der Auflösung Ihres geografischen sensor und die Distanz-Metrik, die Sie verwenden - ich würde vorschlagen, die haversine Entfernung für geographische Daten.

Alle, die ich auf dieser Maschine ist ein Alter compiler, also machte ich eine ASCII-version von diesem. Es "zieht" (in ASCII) eine Karte - Punkte sind Punkte, X ist, wo die wirkliche Quelle ist, G ist, wo die vermutete Quelle ist. Wenn die beiden sich überschneiden, nur X angezeigt.

Beispiele (SCHWIERIGKEIT 1,5 und 3 entsprechend):

Die Punkte erzeugt werden, durch die Auswahl von einem zufälligen Punkt als Quelle, dann nach dem Zufallsprinzip verteilen Punkte, so dass Sie eher näher an der Quelle.

DIFFICULTY ist eine floating-point-Konstanten regelt, dass der Ausgangspunkt generation - wie viel eher die Punkte, um näher an die Quelle - wenn es 1 oder weniger, das Programm sollte in der Lage sein, zu erraten, die genaue Quelle, oder ganz in der Nähe. Bei 2,5, sollte es immer noch ziemlich anständig. Bei 4+, wird es beginnen, zu erraten, noch schlimmer, aber ich denke, dass es immer noch Vermutungen besser als ein Mensch.

Könnte es sein, optimiert durch die Verwendung von binären Suche über die X -, dann Y - dies würde die Vermutung schlimmer, wäre aber viel, viel schneller. Oder beginnen Sie mit größeren Blöcken, dann die Spaltung der beste block weiter (oder der beste block und die 8 umgebenden). Für eine höhere Auflösung-system, für die diese erforderlich wären. Dies ist ein ziemlich naiver Ansatz, obwohl, aber es scheint zu funktionieren gut in einem 80x24-system. 😀

#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);

} 

Einer einfachen Mischung aus Modell scheint zu funktionieren ziemlich gut für dieses problem.

Im Allgemeinen, um einen Punkt minimiert, dass der Abstand zu allen anderen Punkten in einem Datensatz, können Sie einfach den Mittelwert. In diesem Fall, werden Sie wollen, um einen Punkt zu finden, minimiert die Distanz von a Teilmenge der Punkte konzentriert. Wenn Sie postulieren, dass ein Punkt kommen entweder aus der konzentrierten Menge der Punkte von Interesse oder aus einer diffusen Satz von hintergrund-Punkte, dann ergibt dies eine Mischung Modell.

Ich habe einige python-code unten. Die konzentrierten Bereich wird modelliert, indem eine hoch-Präzisions-normal-Verteilung und der hintergrund Punkt modelliert werden, indem entweder eine niedrige Genauigkeit Normalverteilung oder eine Gleichverteilung über eine bounding-box, die auf den Datensatz (es ist eine code-Zeile, die auskommentiert werden können, um wechseln Sie zwischen diesen Optionen). Auch, Mischung Modelle können etwas instabil, so läuft der EM-Algorithmus ein paar mal mit zufälligen Anfangsbedingungen und die Wahl der Lauf mit dem höchsten log-likelihood bessere Ergebnisse liefert.

Wenn Sie eigentlich auf der Suche bei Flugzeugen, dann hinzufügen irgendeine Art von zeitabhängigen Dynamik wird wahrscheinlich verbessern Sie Ihre Fähigkeit zu schließen, die home base immens.

Ich würde auch vorsichtig sein, Rossimo Formel, denn es enthält einige ziemlich starke Annahmen über Kriminalität-Distributionen.

#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()

Zuerst möchte ich sagen, dass ich mit Vorliebe deine Methode illustrieren und erklären das problem ..

Wenn ich in deinen Schuhen wäre, würde ich gehen für einen Dichte-basierter Algorithmus wie DBSCAN
und dann nach der Clusterbildung werden die Gebiete und die Beseitigung der Lärm Punkte wenige Bereiche (Auswahl) bleiben .. dann werde ich den cluster mit der höchsten Dichte der Punkte und Berechnung der durchschnittlichen Punkt und finden Sie die nächste wirkliche Punkt zu . getan, den Platz gefunden! :).

Grüße,

Summe über die Distanzen. Nehmen Sie den Punkt mit der kleinsten summiert Entfernung.

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