Los datos tienen dos tendencias; ¿Cómo extraer líneas de tendencia independientes?

Tengo un conjunto de datos que no está ordenado de ninguna manera en particular, pero cuando se traza claramente tiene dos tendencias distintas. Una regresión lineal simple realmente no sería adecuada aquí debido a la clara distinción entre las dos series. ¿Hay una manera simple de obtener las dos líneas de tendencia lineales independientes?

Para el registro, estoy usando Python y estoy bastante cómodo con la programación y el análisis de datos, incluido el aprendizaje automático, pero estoy dispuesto a saltar a R si es absolutamente necesario.

Para resolver su problema, un buen enfoque es definir un modelo probabilístico que coincida con los supuestos sobre su conjunto de datos. En su caso, probablemente desee una mezcla de modelos de regresión lineal. Puede crear un modelo de "mezcla de regresores" similar a un modelo de mezcla gaussiana asociando diferentes puntos de datos con diferentes componentes de la mezcla.

He incluido un código para que comiences. El código implementa un algoritmo EM para una mezcla de dos regresores (debería ser relativamente fácil de extender a mezclas más grandes). El código parece ser bastante robusto para conjuntos de datos aleatorios. Sin embargo, a diferencia de la regresión lineal, los modelos de mezcla tienen objetivos no convexos, por lo que para un conjunto de datos real, es posible que deba ejecutar algunas pruebas con diferentes puntos de partida aleatorios.

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

En otra parte de este hilo, user1149913 proporciona excelentes consejos (define un modelo probabilístico) y codifica un enfoque poderoso (estimación EM). Quedan dos cuestiones por abordar:

1. Cómo hacer frente a las desviaciones del modelo probabilístico (que son muy evidentes en los datos de 2011-2012 y algo evidentes en las ondulaciones de los puntos menos inclinados).

2. Cómo identificar buenos valores iniciales para el algoritmo EM (o cualquier otro algoritmo).

Para abordar el n. ° 2, considere usar una transformación Hough . Este es un algoritmo de detección de características que, para encontrar extensiones lineales de características, se puede calcular de manera eficiente como una transformación de Radón .

Conceptualmente, la transformación de Hough representa conjuntos de líneas. Una línea en el plano se puede parametrizar por su pendiente, , y su distancia, , desde un origen fijo. Un punto en este sistema de coordenadas tanto designa una sola línea. Cada punto en el diagrama 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 las entidades en el diagrama original caen a lo largo de una línea común, o lo suficientemente cerca de una, entonces las colecciones de curvas que producen en la transformación de Hough tienden a tener una intersección común correspondiente a esa línea común. Al encontrar estos puntos de mayor intensidad en la transformación de Hough, podemos leer buenas soluciones al problema original.y x , $x$$y$$x,y$

Para comenzar con estos datos, primero recorté los elementos auxiliares (ejes, marcas y etiquetas) y, por si acaso, recorté los puntos obviamente externos en la parte inferior derecha y los rocié a lo largo del eje inferior. (Cuando ese material no se recorta, el procedimiento aún funciona bien, ¡pero también detecta los ejes, los marcos, las secuencias lineales de ticks, las secuencias lineales de etiquetas e incluso los puntos que se encuentran esporádicamente en el eje inferior!)

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 están en Mathematica ).

A cada punto de esta imagen corresponde 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 se manifieste visualmente el sentido en que la pregunta es un problema de agrupación de líneas : la transformación de Hough lo reduce a un problema de agrupación de puntos , al que podemos aplicar cualquier método de agrupación que nos guste.

En este caso, la agrupación es tan clara que basta el simple procesamiento posterior de la transformación de Hough. Para identificar ubicaciones de mayor intensidad en la transformación, aumenté el contraste y difuminé la transformación en un radio de aproximadamente 1%: eso es comparable a los diámetros de los puntos de la trama en la imagen original.

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

Limitar el resultado lo redujo a dos pequeñas burbujas cuyos centroides identifican razonablemente los puntos de mayor intensidad: estos estiman las líneas ajustadas.

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

(El umbral de se encontró empíricamente: produce solo dos regiones y la más pequeña de las dos es casi lo más pequeña posible).$0.777$

El lado izquierdo de la imagen corresponde a una dirección de 0 grados (horizontal) y, al mirar de izquierda a derecha, ese ángulo aumenta linealmente a 180 grados. Interpolando, calculo que los dos blobs están centrados a 19 y 57.1 grados, respectivamente. También podemos leer las intersecciones de las posiciones verticales de los blobs. Esta información produce los ajustes iniciales:

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}

De manera similar, uno puede calcular las intersecciones correspondientes a estas pendientes, dando estos ajustes:

(La línea roja corresponde al pequeño punto rosa en la imagen anterior y la línea azul corresponde a la gota de aguamarina más grande).

En gran medida, este enfoque ha abordado automáticamente el primer problema: las desviaciones de la linealidad difuminan los puntos de mayor intensidad, pero generalmente no los cambian demasiado. Los puntos periféricos francos contribuirán con el ruido de bajo nivel a lo largo de la transformación de Hough, que desaparecerá durante los procedimientos de postprocesamiento.

En este punto, se pueden proporcionar estas estimaciones como valores iniciales para el algoritmo EM o para un minimizador de probabilidad (que, dadas buenas estimaciones, convergerán rápidamente). Sin embargo, sería mejor utilizar un estimador de regresión robusto, como los mínimos cuadrados reponderados de forma iterativa . Es capaz de proporcionar un peso de regresión a cada punto. Los pesos bajos indican que un punto no "pertenece" a una línea. Explote estos pesos, si lo desea, para asignar cada punto a su línea adecuada. Luego, una vez clasificados los puntos, puede usar mínimos cuadrados ordinarios (o cualquier otro procedimiento de regresión) por separado en los dos grupos de puntos.

Encontré esta pregunta vinculada a otra pregunta . Realmente hice investigación académica sobre este tipo de problema. Por favor, compruebe mi respuesta "¿Raíz cuadrada mínima"? Un método de ajuste con mínimos mínimos para más detalles.

El enfoque basado en la transformación Hough de Whuber es una muy buena solución para escenarios simples como el que usted proporcionó. Trabajé en escenarios con datos más complejos, como este:

Mis coautores y yo denotamos esto como un problema de "asociación de datos". Cuando intenta resolverlo, el problema principal suele ser combinatorio debido a la cantidad exponencial de posibles combinaciones de datos.

Tenemos una publicación " Mezclas superpuestas de procesos gaussianos para el problema de asociación de datos " donde abordamos el problema general de las curvas N con una técnica iterativa, dando muy buenos resultados. Puede encontrar el código de Matlab vinculado en el documento.

[Actualización] Una implementación de Python de la técnica OMGP se puede encontrar en la biblioteca GPClust .

Tengo otro documento en el que relajamos el problema para obtener un problema de optimización convexo, pero aún no se ha aceptado su publicación. Es específico para 2 curvas, por lo que funcionaría perfectamente en sus datos. Déjame saber si estás interesado.

user1149913 tiene una excelente respuesta (+1), pero me parece que su recopilación de datos se desmoronó a fines de 2011, por lo que tendría que cortar esa parte de sus datos y luego ejecutar las cosas varias veces con diferentes aleatorias coeficientes iniciales para ver lo que obtienes.

Una manera sencilla de hacer las cosas sería separar los datos en dos conjuntos a simple vista, luego usar cualquier técnica de modelo lineal a la que esté acostumbrado. En R, sería la lmfunción.

O ajuste dos líneas a simple vista. En R usarías ablinepara hacer esto.

Los datos están confusos, tienen valores atípicos y se desmoronan al final, pero a simple vista tiene dos líneas bastante obvias, por lo que no estoy seguro de que un método elegante valga la pena.