De datos tiene dos tendencias; cómo extraer independiente de las líneas de tendencia?

Question

Tengo un conjunto de datos que no está ordenado de una manera particular, pero cuando se trazan claramente tiene dos claras tendencias. Una regresión lineal simple no serían realmente adecuado aquí, debido a la clara distinción entre las dos series. Hay una manera sencilla de que los dos independientes lineal de las líneas de tendencia?

Para el registro estoy usando Python y me siento bastante cómodo con la programación y el análisis de datos , incluyendo el aprendizaje de máquina, pero estoy dispuesto a saltar sobre R si es absolutamente necesario.

Answer 1

Para solucionar su problema, un buen enfoque es definir un modelo probabilístico que coincide con la hipótesis acerca de su conjunto de datos. En su caso, usted probablemente querrá una mezcla de los modelos de regresión lineal. Usted puede crear una "mezcla de regresores" modelo similar a un modelo de mezcla de gaussianas por la asociación de puntos de datos diferentes con diferentes componentes de la mezcla.

He incluido algunas código para empezar. El código que implementa un algoritmo EM para una mezcla de dos regresores (que debería ser relativamente fácil ampliar a mayor mezclas). El código parece ser bastante robusto para conjuntos de datos aleatorios. Sin embargo, a diferencia de la regresión lineal, modelos de mezcla no-convexas de los objetivos, por lo que para un verdadero conjunto de datos, puede que tenga que ejecutar un par de ensayos aleatorios diferentes puntos de partida.

import numpy as np
import matplotlib.pyplot as plt 
import scipy.linalg as lin

#generate some random data
N=100
x=np.random.rand(N,2)
x[:,1]=1

w=np.random.rand(2,2)
y=np.zeros(N)

n=int(np.random.rand()*N)
y[:n]=np.dot(x[:n,:],w[0,:])+np.random.normal(size=n)*.01
y[n:]=np.dot(x[n:,:],w[1,:])+np.random.normal(size=N-n)*.01


rx=np.ones( (100,2) )
r=np.arange(0,1,.01)
rx[:,0]=r

#plot the random dataset
plt.plot(x[:,0],y,'.b')
plt.plot(r,np.dot(rx,w[0,:]),':k',linewidth=2)
plt.plot(r,np.dot(rx,w[1,:]),':k',linewidth=2)

# regularization parameter for the regression weights
lam=.01

def em():
    # mixture weights
    rpi=np.zeros( (2) )+.5

    # expected mixture weights for each data point
    pi=np.zeros( (len(x),2) )+.5

    #the regression weights
    w1=np.random.rand(2)
    w2=np.random.rand(2)

    #precision term for the probability of the data under the regression function 
    eta=100

    for _ in xrange(100):
        if 0:
            plt.plot(r,np.dot(rx,w1),'-r',alpha=.5)
            plt.plot(r,np.dot(rx,w2),'-g',alpha=.5)

        #compute lhood for each data point
        err1=y-np.dot(x,w1)
        err2=y-np.dot(x,w2)
        prbs=np.zeros( (len(y),2) )
        prbs[:,0]=-.5*eta*err1**2
        prbs[:,1]=-.5*eta*err2**2

        #compute expected mixture weights
        pi=np.tile(rpi,(len(x),1))*np.exp(prbs)
        pi/=np.tile(np.sum(pi,1),(2,1)).T

        #max with respect to the mixture probabilities
        rpi=np.sum(pi,0)
        rpi/=np.sum(rpi)

        #max with respect to the regression weights
        pi1x=np.tile(pi[:,0],(2,1)).T*x
        xp1=np.dot(pi1x.T,x)+np.eye(2)*lam/eta
        yp1=np.dot(pi1x.T,y)
        w1=lin.solve(xp1,yp1)

        pi2x=np.tile(pi[:,1],(2,1)).T*x
        xp2=np.dot(pi2x.T,x)+np.eye(2)*lam/eta
        yp2=np.dot(pi[:,1]*y,x)
        w2=lin.solve(xp2,yp2)

        #max wrt the precision term
        eta=np.sum(pi)/np.sum(-prbs/eta*pi)

        #objective function - unstable as the pi's become concentrated on a single component
        obj=np.sum(prbs*pi)-np.sum(pi[pi>1e-50]*np.log(pi[pi>1e-50]))+np.sum(pi*np.log(np.tile(rpi,(len(x),1))))+np.log(eta)*np.sum(pi)
        print obj,eta,rpi,w1,w2

        try:
            if np.isnan(obj): break
            if np.abs(obj-oldobj)<1e-2: break
        except:
            pass

        oldobj=obj

    return w1,w2


#run the em algorithm and plot the solution
rw1,rw2=em()
plt.plot(r,np.dot(rx,rw1),'-r')
plt.plot(r,np.dot(rx,rw2),'-g')

plt.show()

Answer 2

En otra parte de este hilo, user1149913 proporciona un gran consejo (definir un modelo probabilístico) y el código para un enfoque de gran alcance (EM estimación). Dos problemas por resolver:

1. Cómo lidiar con salidas desde el modelo probabilístico (que son muy evidentes en el 2011-2012 datos y algo evidente en las ondulaciones de las menos inclinadas puntos).

2. Cómo identificar los valores de partida para el algoritmo EM (o cualquier otro algoritmo).

A la dirección #2, considere el uso de una transformación de Hough. Esta es una característica algoritmo de detección que, para encontrar lineal tramos de características, de manera eficiente puede ser calculada como el Radón transformar.

Conceptualmente, la transformación de Hough representa conjuntos de líneas. Una línea en el plano puede ser parametrizado por su pendiente, $x$, y su distancia, $y$, desde un punto fijo de origen. Un punto en este $x,$ y sistema de coordenadas de lo que designa una sola línea. Cada punto en la trama original determina un lápiz de líneas que pasan por ese punto: este lápiz aparece como una curva en la transformación de Hough. Cuando los elementos en la trama original caída a lo largo de una línea común, o lo suficientemente cerca como para uno, entonces las colecciones de curvas que producen en la transformación de Hough tienden a tener un común de intersección correspondiente a la línea común. Al encontrar estos puntos de mayor intensidad en la transformación de Hough, podemos leer buenas soluciones para el problema original.

Para empezar con estos datos,la primera vez me quede afuera de la auxiliar de cosas (ejes, las marcas de graduación y etiquetas) y en buena medida recortado el obviamente periféricas de puntos en la parte inferior derecha y la esparcieron a lo largo de la parte inferior del eje. (Cuando eso no es recortado, el procedimiento todavía funciona bien, sino que también detecta los ejes, los marcos, las secuencias lineales de las garrapatas, las secuencias lineales de etiquetas, e incluso los puntos de mentir de forma esporádica en la parte inferior del eje!)

img = Import["http://i.stack.imgur.com/SkEm3.png"]
i = ColorNegate[Binarize[img]]
crop2 = ImageCrop[ImageCrop[i, {694, 531}, {Left, Bottom}], {565, 467}, {Right, Top}]

(Este y el resto del código en Mathematica.)

Para cada punto en esta imagen corresponde a un rango estrecho de curvas en la transformación de Hough, visible aquí. Son ondas sinusoidales:

hough2 = Radon[crop2, Method -> "Hough"]  // ImageAdjust

Esto hace que visualmente manifiesto el sentido en el que la pregunta es una línea de la agrupación problema: la transformación de Hough, se reduce a un punto de la agrupación problema, para el cual se puede aplicar cualquier método de agrupación que nos gusta.

En este caso, la agrupación está tan claro que el post-procesamiento simple de la transformación de Hough bastado. Identificar las zonas de mayor intensidad de la transformación, he aumentado el contraste y borrosa de la transformación en un radio de alrededor de 1%: eso es comparable a la de los diámetros de los puntos de la trama en la imagen original.

blur = ImageAdjust[Blur[ImageAdjust[hough2, {1, 0}], 8]]

Umbralización el resultado condujo a dos pequeñas manchas cuyos centroides razonablemente identificar los puntos de mayor intensidad: estos estimar el conjunto de líneas.

comp = MorphologicalComponents[blur, 0.777]) // Colorize

(El umbral de $0.777$ se ha encontrado empíricamente: se produce sólo dos regiones y el menor de los dos es casi tan pequeño como sea posible.)

El lado izquierdo de la imagen corresponde a una dirección de 0 grados (horizontal) y, como vemos, de izquierda a derecha, que el ángulo aumenta linealmente a 180 grados. La interpolación, calculo que las dos manchas se centran en 19 y 57.1 grados, respectivamente. También podemos leer en las intersecciones de las posiciones verticales de los blobs. Esta información se obtiene información inicial se ajusta:

width = ImageDimensions[blur][[1]];
slopes =  Module[{x, y, z}, ComponentMeasurements[comp, "Centroid"] /. 
          Rule[x_, {y_, z_}] :>  Round[((y - 1/2)/(width - 1))  180., 0.1]
  ]

{19., 57.1}

En una manera similar, uno puede calcular las intersecciones correspondientes a estas laderas, dando a estos ajustes:

(La línea roja corresponde a el pequeño punto de color rosa en la foto anterior y la línea azul corresponde a la mayor aqua blob.)

En gran medida, este enfoque tiene repartidas automáticamente con el primer problema: desviaciones de la linealidad citología por los puntos de mayor intensidad, pero normalmente no cambio mucho. Francamente periféricas puntos contribuirá bajo nivel de ruido durante la transformación de Hough, que desaparecerán durante el post-procesamiento de los procedimientos.

En este punto, se puede proporcionar estas estimaciones como valores de partida para el algoritmo EM o para una probabilidad minimizer (que, dado el buen estimaciones, convergerá rápidamente). Mejor, sin embargo, sería el uso de un estimador de regresión robusta como de forma iterativa reponderadas de los mínimos cuadrados. Es capaz de proporcionar una regresión de peso para cada punto. Pesos bajos indican un punto de no "pertenecen" a una línea. La explotación de estos pesos, si se desea, para asignar a cada punto de su línea. Luego, después de haber clasificado los puntos, puede utilizar mínimos cuadrados ordinarios (o de cualquier otro procedimiento de regresión) por separado en los dos grupos de puntos.

Answer 3

He encontrado esta pregunta vinculada a otra pregunta. De hecho, me hizo la investigación académica en este tipo de problema. Por favor revise mi respuesta "Menos raíz cuadrada de" ajuste? Un método de ajuste con múltiples mínimos para obtener más detalles.

whuber del Hough transformar el enfoque basado en una muy buena solución para escenarios simples como el que usted le dio. He trabajado en escenarios con datos más complejos, como este:

Mis coautores y yo denota una "asociación de datos" problema. Cuando intenta resolverlo, el principal problema suele ser la combinatoria debido a la exponencial de la cantidad de posibles combinaciones de datos.

Tenemos una publicación "la Superposición de Mezclas de Gaussianas Procesos para los datos de la asociación de problema", donde abordamos el problema general de la N curvas con una técnica iterativa, dando muy buenos resultados. Usted puede encontrar el código de Matlab vinculados en el papel (lo siento, no R ni Python).

Tengo otro papel donde nos relajamos el problema de obtener un problema de optimización convexa, pero no ha sido aceptado para su publicación todavía. Es específico para 2 curvas, por lo que funciona a la perfección en sus datos. Déjeme saber si usted está interesado.

Answer 4

user1149913 tiene una excelente respuesta (+1), pero a mi me parece que la recolección de datos se desplomó a finales de 2011, así que vas a tener que cortar la parte de los datos y, a continuación, ejecutar las cosas un par de veces con diferentes al azar a partir de los coeficientes para ver lo que se obtiene.

Una sencilla manera de hacer las cosas, sería la de separar los datos en dos conjuntos por el ojo, a continuación, utilice cualquiera que sea el modelo lineal de la técnica a la que estás acostumbrado. En R, sería el lm función.

O colocar dos líneas por los ojos. En R utilizaría abline a ello.

Los datos están mezcladas, tiene valores atípicos, y se cae en la final, sin embargo, por-ojo tiene dos bastante obvio líneas, por lo que no estoy seguro de que una fantasía método es la pena.