algorithm - para - ¿Cómo puedo encontrar el centro de un grupo de puntos de datos?
hay grupos en instagram (14)
En astrofísica usamos el llamado "radio de media masa". Dada una distribución y su centro, el radio de media masa es el radio mínimo de un círculo que contiene la mitad de los puntos de su distribución.
Esta cantidad es una longitud característica de una distribución de puntos.
¡Si quieres que la casa del helicóptero sea donde los puntos están concentrados al máximo, es el punto que tiene el radio mínimo de media masa!
Mi algoritmo es el siguiente: para cada punto, calcula este radio de media masa centrando la distribución en el punto actual. El "hogar" del helicóptero será el punto con el radio mínimo de la mitad de la masa.
Lo he implementado y el centro computado es 42.149994 -88.133698 (que está en Chicago). También he usado el 0.2 de la masa total en vez del 0.5 (la mitad) generalmente usado en Astrofísica.
Este es mi algoritmo (en python) que encuentra el hogar del helicóptero:
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]
Supongamos que he trazado la posición de un helicóptero todos los días durante el año pasado y se me ocurrió el siguiente mapa:
Cualquier persona que vea esto podría decirme que este helicóptero tiene su base en Chicago.
¿Cómo puedo encontrar el mismo resultado en el código?
Estoy buscando algo como esto:
$geoCodeArray = array([GET=http://pastebin.com/grVsbgL9]);
function findHome($geoCodeArray) {
// magic
return $geoCode;
}
En última instancia generando algo como esto:
ACTUALIZACIÓN: Dataset de muestra
Aquí hay un mapa con un conjunto de datos de muestra: http://batchgeo.com/map/c3676fe29985f00e1605cd4f86920179
Aquí hay un pastebin de 150 geocódigos: http://pastebin.com/grVsbgL9
Lo anterior contiene 150 geocódigos. Los primeros 50 están en algunos grupos cerca de Chicago. El resto está disperso por todo el país, incluidos algunos pequeños grupos en Nueva York, Los Ángeles y San Francisco.
Tengo alrededor de un millón de (seriamente) conjuntos de datos como este que necesitaré para iterar e identificar el "hogar" más probable. Tu ayuda es muy apreciada.
ACTUALIZACIÓN 2: Avión cambiado a helicóptero
El concepto de avión atraía demasiada atención hacia los aeropuertos físicos. Las coordenadas pueden estar en cualquier parte del mundo, no solo en los aeropuertos. Supongamos que es un súper helicóptero no limitado por la física, el combustible o cualquier otra cosa. Puede aterrizar donde quiera. ;)
Encuentra el punto con el mayor estimado de densidad .
Debería ser bastante sencillo. Use un radio de núcleo que cubra aproximadamente un gran aeropuerto de diámetro. Un kernel 2D gaussiano o Epanechnikov debería estar bien.
http://en.wikipedia.org/wiki/Multivariate_kernel_density_estimation
Esto es similar a calcular un Mapa de Heap: http://en.wikipedia.org/wiki/Heat_map y luego encontrar el lugar más brillante allí. Excepto que calcula el brillo de inmediato.
Por diversión, leí una muestra del 1% de Geocoordinadas de DBpedia (es decir, Wikipedia) en ELKI, la proyecté en el espacio 3D y habilité la superposición de estimación de densidad (oculta en el menú de diagrama de dispersión de los visualizadores). Puedes ver que hay un punto de acceso en Europa y, en menor medida, en los EE. UU. El punto de acceso en Europa es Polonia, creo. La última vez que lo comprobé, aparentemente alguien había creado un artículo de Wikipedia con geocoordinados para casi cualquier ciudad de Polonia. Desafortunadamente, el visualizador ELKI no le permite acercarse, girar o reducir el ancho de banda del kernel para encontrar visualmente el punto más denso. Pero es sencillo implementarlo usted mismo; Es probable que tampoco necesite ir al espacio 3D, sino que solo puede usar latitudes y longitudes.
La estimación de la densidad del núcleo debe estar disponible en toneladas de aplicaciones. El que está en R es probablemente mucho más poderoso. Hace poco descubrí este mapa de calor en ELKI, así que supe cómo acceder a él rápidamente. Consulte, por ejemplo, http://stat.ethz.ch/R-manual/R-devel/library/stats/html/density.html para ver una función R relacionada.
En tus datos, en R, prueba por ejemplo:
library(kernSmooth)
smoothScatter(data, nbin=512, bandwidth=c(.25,.25))
Esto debería mostrar una fuerte preferencia por 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]])
cede [1] 42.14697 -88.09508 , que está a menos de 10 millas del aeropuerto de Chicago.
Para obtener mejores coordenadas prueba:
- volver a ejecutar en un área de 20x20 millas alrededor de las coordenadas estimadas
- un KDE sin contenedor en esa área
- mejor selección de ancho de banda con
dpik - mayor resolución de cuadrícula
Esto se puede resolver encontrando una superficie de peligro . Ver la Fórmula de Rossmo .
Este es el problema de los depredadores . Dado un conjunto de cadáveres ubicados geográficamente, ¿dónde está la guarida del depredador? La fórmula de Rossmo resuelve este problema.
La siguiente solución funciona incluso si los puntos están dispersos por toda la Tierra, al convertir la latitud y la longitud en coordenadas cartesianas. Hace un tipo de KDE (estimación de la densidad del kernel), pero en un primer paso la suma de los kernels se evalúa solo en los puntos de datos. El núcleo debe ser elegido para adaptarse al problema. En el siguiente código, es lo que podría llamar en broma / presuntuosamente un Trossian, es decir, 2-d² / h² para d≤h y h² / d² para d> h (donde d es la distancia euclidiana yh es el "ancho de banda" $global_kernel_radius ), pero también podría ser un Gaussian (e -d² / 2h² ), un kernel Epanechnikov (1-d² / h² para d <h, 0 de lo contrario), u otro kernel. Un segundo pase opcional refina la búsqueda localmente, ya sea sumando un kernel independiente en una grilla local, o calculando el centroide, en ambos casos en un entorno definido por $local_grid_radius .
En esencia , cada punto suma todos los puntos que tiene alrededor (incluido él mismo), ponderándolos más si están más cerca (por la curva de campana), y también $w_arr matriz de peso opcional $w_arr . El ganador es el punto con la suma máxima. Una vez que se ha encontrado al ganador, el "hogar" que estamos buscando se puede encontrar repitiendo el mismo proceso localmente alrededor del ganador (usando otra curva de campana), o puede estimarse como el "centro de masa" de todos los puntos dentro de un radio dado del ganador, donde el radio puede ser cero.
El algoritmo debe adaptarse al problema seleccionando los núcleos adecuados, eligiendo cómo refinar la búsqueda localmente y ajustando los parámetros. Para el conjunto de datos de ejemplo, el kernel Trossian para la primera pasada y el kernel Epanechnikov para la segunda pasada, con los 3 radios configurados en 30 mi y un paso de 1 mil en la cuadrícula podrían ser un buen punto de inicio, pero solo si los dos subgrupos Los grupos de Chicago deben ser vistos como un grupo grande. De lo contrario se deben elegir radios más pequeños.
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);
}
El hecho de que las distancias sean euclidianas y no de gran círculo debería tener efectos insignificantes para la tarea en cuestión. Calcular las distancias del círculo máximo sería mucho más engorroso y haría que el peso de los puntos muy lejanos sea significativamente menor, pero estos puntos ya tienen un peso muy bajo. En principio, un kernel diferente podría lograr el mismo efecto. Los núcleos que tienen un punto de corte completo más allá de cierta distancia, como el núcleo Epanechnikov, no tienen este problema en absoluto (en la práctica).
La conversión entre lat, lng y x, y, z para el dato WGS84 se da exactamente (aunque sin garantía de estabilidad numérica) más como referencia que por una necesidad real. Si se debe tener en cuenta la altura o si se necesita una conversión de retroceso más rápida, consulte el artículo de Wikipedia .
El kernel de Epanechnikov, además de ser "más local" que los núcleos de Gauss y Trossian, tiene la ventaja de ser el más rápido para el segundo ciclo, que es O (ng), donde g es el número de puntos de la cuadrícula local y puede También se empleará en el primer bucle, que es O (n²), si n es grande.
Puede adaptar fácilmente la fórmula de Rossmo, citada por Tyler Durden a su caso, con algunas notas simples:
La formula :
Esta fórmula proporciona algo parecido a la probabilidad de presencia de la operación base para un depredador o un asesino en serie. En tu caso, podría dar la probabilidad de que una base esté en un cierto punto. Más adelante explicaré cómo usarlo. U puede escribirlo de esta manera:
Proba (base en el punto A) = Suma {en todos los puntos} (Phi / (dist ^ f) + (1-Phi) (B * (gf)) / (2B-dist) ^ g)
Usando la distancia euclidiana
Desea una distancia euclidiana y no la de Manhattan porque un avión o un helicóptero no está obligado a la carretera / calles. Por lo tanto, usar la distancia euclidiana es la forma correcta, si está siguiendo un avión y no un asesino en serie. Entonces, "dist" en la fórmula es la distancia euclidiana entre el punto de prueba y el punto considerado
Tomando la variable razonable B
La variable B se utilizó para representar la regla "un asesino razonablemente inteligente no matará a su vecino". En su caso, el testamento también se aplicó porque nadie usa un avión / roflcopter para llegar a la siguiente esquina. Podemos suponer que el viaje mínimo es, por ejemplo, de 10 km o algo razonable cuando se aplica a su caso.
Factor exponencial f
El factor f se usa para agregar un peso a la distancia. Por ejemplo, si todas las manchas se encuentran en un área pequeña, podría desear un gran factor f porque la probabilidad de que el aeropuerto / base / HQ disminuya rápidamente si todos sus puntos de datos están en el mismo sector. g funciona de manera similar, permite elegir el tamaño de la "base no es probable que esté justo al lado del lugar"
Factor phi:
Nuevamente, este factor debe determinarse utilizando su conocimiento del problema. Permite elegir el factor más preciso entre "la base está cerca de los puntos" y "no usaré el plano para hacer 5 m" si, por ejemplo, cree que la segunda es casi irrelevante, puede configurar Phi en 0.95 (0<Phi<1)Si ambos son interesantes phi estarán alrededor de 0.5
Cómo implementarlo como algo útil:
Primero, debe dividir su mapa en cuadraditos: entrelazar el mapa (como invisal) (cuanto más pequeños sean los cuadrados, más preciso será el resultado (en general)) y luego usar la fórmula para encontrar la ubicación más probable. De hecho, la malla es solo una matriz con todas las ubicaciones posibles. (Si quieres ser preciso, aumentas la cantidad de puntos posibles, pero requerirá más tiempo de computación y PhP no es muy conocido por su increíble velocidad)
Algoritmo
//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''s 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))
Espero que te ayude
Puede usar DBSCAN para esa tarea.
DBSCAN es un clúster basado en densidad con una noción de ruido. Necesitas dos parámetros:
Primero, la cantidad de puntos que un grupo debe tener como mínimo "minpoints" mínimos "minpoints" . Y segundo, un parámetro de barrio llamado "epsilon" que establece un umbral de distancia a los puntos circundantes que deben incluirse en el clúster.
Todo el algoritmo funciona así:
- Comience con un punto arbitrario en su conjunto que aún no se haya visitado
- Recuperar puntos del vecindario epsilon marcar todos como visitados
- Si ha encontrado suficientes puntos en este vecindario (> parámetro de minpoints), inicia un nuevo grupo y asigna esos puntos. Ahora recurse en el paso 2 nuevamente para cada punto de este clúster.
- si no tiene, declare este punto como ruido
- vuelve a repetir hasta que hayas visitado todos los puntos
Es realmente simple de implementar y hay muchos frameworks que soportan este algoritmo. Para encontrar la media de su grupo, simplemente puede tomar la media de todos los puntos asignados de su vecindario.
Sin embargo, a diferencia del método que @TylerDurden propone, esto necesita una parametrización, por lo que debe encontrar algunos parámetros ajustados a mano que se ajusten a su problema.
En su caso, puede tratar de establecer los minpoints al 10% de su total de puntos si es probable que el avión permanezca el 10% del tiempo que rastrea en un aeropuerto. El parámetro de densidad epsilon depende de la resolución de su sensor geográfico y de la métrica de distancia que utiliza. Yo sugeriría la distancia más grande para los datos geográficos.
Todo lo que tengo en esta máquina es un viejo compilador, así que hice una versión ASCII de esto. "Dibuja" (en ASCII) un mapa: los puntos son puntos, X es donde está la fuente real, G es donde está la fuente adivinada. Si los dos se superponen, solo se muestra X
Ejemplos (DIFICULTAD 1.5 y 3 respectivamente):
Los puntos se generan seleccionando un punto aleatorio como fuente, luego distribuyendo puntos aleatoriamente, lo que hace que sea más probable que estén más cerca de la fuente.
DIFFICULTY es una constante de punto flotante que regula la generación de puntos inicial (cuanto más probable es que los puntos estén más cerca de la fuente), si es 1 o menos, el programa debe poder adivinar la fuente exacta o muy cerca. En 2.5, todavía debería ser bastante decente. Con 4+, comenzará a adivinar peor, pero creo que todavía adivina mejor de lo que lo haría un humano.
Podría optimizarse usando la búsqueda binaria sobre X, luego Y; esto empeoraría la suposición, pero sería mucho, mucho más rápido. O comenzando con bloques más grandes, luego dividiendo mejor el mejor bloque (o el mejor bloque y los 8 que lo rodean). Para un sistema de mayor resolución, uno de estos sería necesario. Sin embargo, este es un enfoque bastante ingenuo, pero parece funcionar bien en un sistema de 80x24. :RE
#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);
}
Un modelo de mezcla simple parece funcionar bastante bien para este problema.
En general, para obtener un punto que minimiza la distancia a todos los otros puntos en un conjunto de datos, puede simplemente tomar la media. En este caso, desea buscar un punto que minimice la distancia desde un subconjunto de puntos concentrados. Si postula que un punto puede provenir del conjunto concentrado de puntos de interés o de un conjunto difuso de puntos de fondo, entonces se obtiene un modelo de mezcla.
He incluido algunos códigos Python a continuación. El área concentrada se modela mediante una distribución normal de alta precisión y el punto de fondo se modela mediante una distribución normal de baja precisión o una distribución uniforme sobre un cuadro delimitador en el conjunto de datos (hay una línea de código que puede comentarse cambiar entre estas opciones). Además, los modelos de mezcla pueden ser algo inestables, por lo que ejecutar el algoritmo EM varias veces con condiciones iniciales aleatorias y elegir la ejecución con la mayor probabilidad de registro da mejores resultados.
Si realmente está mirando aviones, entonces agregar algún tipo de dinámica dependiente del tiempo probablemente mejorará su capacidad de inferir inmensamente la base de operaciones.
También desconfiaría de la fórmula de Rossimo porque incluye algunas suposiciones bastante fuertes sobre la distribución de delitos.
#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()
Virtual Earth tiene una muy buena explicación de cómo puedes hacerlo relativamente rápido. También han proporcionado ejemplos de código. Por favor, eche un vistazo a http://soulsolutions.com.au/Articles/ClusteringVirtualEarthPart1.aspx
¿Por qué no algo como esto?
- Para cada punto, calcula la distancia desde todos los demás puntos y suma el total.
- El punto con la suma más pequeña es tu centro.
Tal vez la suma no es la mejor métrica para usar. Posiblemente el punto con más "pequeñas distancias"?
Primero me gustaría expresar mi cariño por su método al ilustrar y explicar el problema ...
Si estuviera en su lugar, buscaría un algoritmo basado en la densidad como DBSCAN y luego, después de agrupar las áreas y eliminar los puntos de ruido , quedarán algunas (opciones) ... luego tomaré el grupo con la mayor densidad de apunta y calcula el punto promedio y encuentra el punto real más cercano . Hecho, encontré el lugar! :).
Saludos,
¿Qué hay de dividir el mapa en muchas zonas y luego encontrar el centro del plano en la zona con más plano? Algoritmo será algo como esto
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
Puede tomar un árbol de expansión mínimo y eliminar los bordes más largos. Los árboles más pequeños te dan el centeroide para buscar. El nombre del algoritmo es k-clustering de enlace único. Hay una publicación aquí: https://stats.stackexchange.com/questions/1475/visualization-software-for-clustering .
Suma sobre las distancias. Toma el punto con la menor distancia sumada.
function () {
for i in points P:
S[i] = 0
for j in points P:
S[i] += distance(P[i], P[j])
return min(S);
}