algorithm - データ ポイントのクラスターの中心を見つけるにはどうすればよいですか?

Question

過去 1 年間、毎日ヘリコプターの位置をプロットし、次のマップを作成したとします。

これを見た人なら誰でも、このヘリコプターがシカゴを拠点としているとわかるでしょう。

コードで同じ結果を見つけるにはどうすればよいですか?

私はこのようなものを探しています:

$geoCodeArray = array([GET=http://pastebin.com/grVsbgL9]);
function findHome($geoCodeArray) {
    // magic
    return $geoCode;
}

最終的に次のようなものを生成します。

更新: サンプル データセット

サンプル データセットを含むマップを次に示します: http://batchgeo.com/map/c3676fe29985f00e1605cd4f86920179

150 個のジオコードのペーストビンは次のとおりです: http://pastebin.com/grVsbgL9

上記には 150 個のジオコードが含まれています。最初の 50 は、シカゴに近いいくつかのクラスターに属しています。残りは、ニューヨーク、ロサンゼルス、サンフランシスコのいくつかの小さなクラスターを含め、全国に散らばっています.

私は、このような (真剣に) 約 100 万のデータセットを持っており、それらを反復処理して、最も可能性の高い「家」を特定する必要があります。よろしくお願いいたします。

更新 2: 飛行機がヘリコプターに切り替わった

飛行機のコンセプトは、物理的な空港に向けられすぎていました。座標は、空港だけでなく、世界中のどこでもかまいません。物理学や燃料などに縛られないスーパー ヘリコプタだとしましょう。好きな場所に着陸できます。;)

Answer 1

次の解決策は、緯度と経度をデカルト座標に変換することにより、ポイントが地球全体に分散している場合でも機能します。これは一種の KDE (カーネル密度推定) を行いますが、最初のパスでは、カーネルの合計はデータ ポイントでのみ評価されます。問題に適合するようにカーネルを選択する必要があります。以下のコードでは、私が冗談めかして/思い切ってトロシアンと呼ぶことができるものです。つまり、d≤h の場合は 2-d²/h²、d>h の場合は h²/d² です (d はユークリッド距離、h は「帯域幅」$global_kernel_radius) 。ですが、ガウス (e ^-d²/2h² )、エパネチニコフ カーネル (d<h の場合は 1-d²/h²、それ以外の場合は 0)、または別のカーネルの場合もあります。オプションの 2 番目のパスでは、ローカル グリッドで独立したカーネルを合計するか、重心を計算することにより、検索をローカルに絞り込みます。$local_grid_radius.

本質的に、各ポイントは、周囲にあるすべてのポイント (それ自体を含む) を合計し、それらが近くにある場合は (ベル カーブによって) 重み付けし、オプションの weight array によって重み付けします$w_arr。勝者は、合計が最大のポイントです。勝者が見つかったら、探している「ホーム」は、同じプロセスを勝者の周りでローカルに繰り返すことで見つけることができます (別のベル曲線を使用)。または、すべてのポイントの「重心」であると推定できます。勝者からの指定された半径内。半径はゼロにすることができます。

適切なカーネルを選択し、ローカルで検索を絞り込む方法を選択し、パラメーターを調整することにより、アルゴリズムを問題に適応させる必要があります。サンプル データセットの場合、1 回目のパスに Trossian カーネル、2 回目のパスに Epanechnikov カーネルを使用し、3 つの半径すべてを 30 マイルに設定し、グリッド ステップを 1 マイルに設定するのが適切な出発点になる可能性があります。シカゴのクラスターは、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 カーネルは、Gaussian および Trossian カーネルよりも「ローカル」であることに加えて、2 番目のループで最速であるという利点があります。これは、g がローカル グリッドのポイントの数である O(ng) であり、次のことができます。 n が大きい場合、O(n²) である最初のループでも使用されます。

Answer 2

これは、危険な表面を見つけることで解決できます。Rossmo の式を参照してください。

これが捕食者の問題です。地理的に配置された一連の死体が与えられた場合、捕食者の隠れ家はどこですか? Rossmo の公式は、この問題を解決します。

Answer 3

推定密度が最大の点を見つけます。

かなり簡単なはずです。直径の大きい空港を大まかにカバーするカーネル半径を使用します。2D Gaussian または Epanechnikov カーネルで問題ありません。

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

これは、ヒープ マップを計算するのと似ています: http://en.wikipedia.org/wiki/Heat_map で、そこで最も明るいスポットを見つけます。ただし、明るさはすぐに計算されます。

楽しみのために、DBpedia (つまりウィキペディア) の Geocoordinates の 1% のサンプルを ELKI に読み込み、それを 3D 空間に投影し、密度推定オーバーレイ (ビジュアライザーの散布図メニューに隠されています) を有効にしました。ヨーロッパにホットスポットがあり、米国にもホットスポットがあることがわかります。ヨーロッパのホットスポットはポーランドだと思います。最後に確認したところ、誰かがポーランドのほぼすべての町の地理座標を使用してウィキペディアの記事を作成していたようです。残念ながら、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
• より高いグリッド解像度

Answer 4

天体物理学では、いわゆる「半質量半径」を使用します。分布とその中心が与えられた場合、半質量半径は、分布の点の半分を含む円の最小半径です。

この量は、点の分布の代表的な長さです。

ヘリコプターのホームをポイントが最大に集中する場所にしたい場合は、半質量半径が最小のポイントになります。

私のアルゴリズムは次のとおりです。各ポイントについて、現在のポイントの分布を中心とするこの半質量半径を計算します。ヘリコプターの「ホーム」は、半質量半径が最小のポイントになります。

私はそれを実装し、計算された中心は42.149994 -88.133698(シカゴにあります) また、天体物理学で通常使用される 0.5(半分) の代わりに、総質量の 0.2 を使用しました。

これは、ヘリコプターのホームを見つける私の (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]

Answer 5

そのタスクにはDBSCANを使用できます。

DBSCAN は、ノイズの概念を使用した密度ベースのクラスタリングです。次の 2 つのパラメーターが必要です。

最初に、クラスターが最低限持つべきポイントの数"minpoints"。"epsilon"次に、クラスターに含める必要がある周囲のポイントまでの距離のしきい値を設定する近傍パラメーターが呼び出されます。

アルゴリズム全体は次のように機能します。

1. まだ訪れていないセット内の任意のポイントから始めます
2. イプシロン近傍からポイントを取得し、すべてを訪問済みとしてマークします
   1. この近傍 (> minpoints パラメーター) で十分なポイントが見つかった場合は、新しいクラスターを開始してそれらのポイントを割り当てます。このクラスター内のすべてのポイントについて、ステップ 2 に戻ります。
   2. 持っていない場合は、この点をノイズとして宣言してください
3. すべてのポイントを訪問するまで、もう一度やり直してください

実装は非常に簡単で、このアルゴリズムをサポートするフレームワークはすでにたくさんあります。クラスターの平均を見つけるには、その近傍から割り当てられたすべてのポイントの平均を取るだけです。

ただし、@TylerDurden が提案する方法とは異なり、これにはパラメーター化が必要です。そのため、問題に合った手動で調整されたパラメーターを見つける必要があります。

あなたの場合、飛行機が空港で追跡する時間の 10% に留まる可能性が高い場合は、最小ポイントを合計ポイントの 10% に設定することができます。密度パラメーターのイプシロンは、地理センサーの解像度と使用する距離メトリックに依存します。地理データにはハーサイン距離をお勧めします。

Answer 6

このマシンには古いコンパイラしかないので、これの ASCII バージョンを作成しました。これは (ASCII で) マップを「描画」します - 点は点、X は実際のソースの場所、G は推測されたソースの場所です。2 つが重複する場合は、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);

} 

Answer 7

マップを多くのゾーンに分割し、最も平面の多いゾーンで平面の中心を見つけてみませんか。アルゴリズムはこのようなものになります

Zones[40]
 foreach Plane in Planes を設定します
   Zones[GetZone(Plane.position)].Add(Plane)

MaxZone.Length() < Zone.Length()の場合、 MaxZone = Zones[0]
 foreach Zone in Zones を
   設定します
           MaxZone = ゾーン

MaxZoneの各平面の中心
を設定します
     中心.X += 平面.X
     中心.Y += 平面.Y
Center.X /= MaxZone.Length
Center.Y /= MaxZone.Length

Answer 8

最初に、問題を図解して説明するあなたの方法への私の好意を表明したいと思います..

私があなたの立場なら、 DBSCANのような密度ベースのアルゴリズムを使用し、領域をクラスタリングしてノイズ ポイントを削除した後 、いくつかの領域 (選択肢) が残ります..次に、密度が最も高いクラスタを取得します。ポイントと平均ポイントを計算し、それに最も近い実際のポイントを見つけます。完了、場所を見つけました！:)。

よろしく、

Answer 9

この問題には、単純な混合モデルがうまく機能するようです。

一般に、データセット内の他のすべてのポイントまでの距離を最小にするポイントを取得するには、平均を取ることができます。この場合、集中したポイントのサブセットからの距離を最小化するポイントを見つける必要があります。ポイントが集中した関心ポイントのセットまたは背景ポイントの拡散セットのいずれかに由来すると仮定すると、混合モデルが得られます。

以下にいくつかのPythonコードを含めました。集中領域は高精度の正規分布によってモデル化され、背景点は低精度の正規分布またはデータセットの境界ボックス上の一様分布のいずれかによってモデル化されます (コメントアウトできるコード行があります)これらのオプションを切り替えます)。また、混合モデルは多少不安定になる可能性があるため、ランダムな初期条件で EM アルゴリズムを数回実行し、対数尤度が最も高い実行を選択すると、より良い結果が得られます。

実際に飛行機を見ているのであれば、ある種の時間依存のダイナミクスを追加すると、ホームベースを推測する能力が大幅に向上する可能性があります。

Rossimo の式には、犯罪分布に関するかなり強力な仮定が含まれているため、私も注意が必要です。

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

Answer 10

最小スパニング ツリーを取得して、最長のエッジを削除できます。小さいツリーは、ルックアップするためのセンターロイドを提供します。アルゴリズム名は、シングルリンク k-クラスタリングです。ここに投稿があります: https://stats.stackexchange.com/questions/1475/visualization-software-for-clustering。

Answer 11

なぜこのようなものではありません：

• 各ポイントについて、他のすべてのポイントからの距離を計算し、合計します。
• 合計が最小の点があなたの中心です。

おそらく、合計は使用するのに最適なメトリックではありません。おそらく、最も「距離が短い」ポイントはありますか?

Answer 12

距離を合計します。距離の合計が最も小さい点を取ります。

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