Les données ont deux tendances; Comment extraire des courbes de tendance indépendantes?

34

J'ai un ensemble de données qui n'est pas ordonné de manière particulière, mais qui présente clairement deux tendances distinctes. Une régression linéaire simple ne conviendrait pas vraiment ici à cause de la distinction claire entre les deux séries. Existe-t-il un moyen simple d’obtenir les deux courbes de tendance linéaires indépendantes?

Pour mémoire, j'utilise Python et je suis assez à l'aise avec la programmation et l'analyse des données, y compris l'apprentissage automatique, mais je suis prêt à passer à R si cela est absolument nécessaire.

entrez la description de l'image ici

jbbiomed
la source
6
La meilleure réponse que j’ai jusqu’à présent est de l’imprimer sur du papier quadrillé et d’utiliser un crayon, une règle et une calculatrice ...
jbbiomed le
Vous pouvez peut-être calculer des pentes par paires et les regrouper en deux "groupes de pentes". Cependant, cela échouera si vous avez deux tendances parallèles.
Thomas Jungblut
1
Je n'ai pas d'expérience personnelle avec cela, mais je pense que les stats modèles mériteraient d'être vérifiés. Statistiquement, une régression linéaire avec une interaction pour groupe serait adéquate (sauf si vous dites que vous avez des données non groupées, auquel cas c'est un peu plus poilu ...)
Matt Parker
1
Malheureusement, il ne s’agit pas de données d’effet, mais de données d’utilisation, et il est clair que l’utilisation de deux systèmes distincts est confondue dans le même ensemble de données. Je veux pouvoir décrire les deux modèles d'utilisation, mais je ne peux pas revenir en arrière et me rappeler des données, car cela représente environ 6 années d'informations collectées par un client.
jbbiomed
2
Juste pour vous assurer: votre client ne dispose pas de données supplémentaires qui indiqueraient quelles mesures proviennent de quelle population? Il s'agit de 100% des données que vous ou votre client possédez ou pouvez trouver. En outre, 2012 semble indiquer que votre collection de données s'est effondrée ou que l'un de vos systèmes ou les deux sont tombés à l'eau. Je me demande si les tendances jusqu'à ce point importent beaucoup.
Wayne

Réponses:

30

Pour résoudre votre problème, une bonne approche consiste à définir un modèle probabiliste qui correspond aux hypothèses sur votre jeu de données. Dans votre cas, vous voulez probablement un mélange de modèles de régression linéaire. Vous pouvez créer un modèle de "mélange de régresseurs" similaire à un modèle de mélange gaussien en associant différents points de données à différentes composantes du mélange.

J'ai inclus du code pour vous aider à démarrer. Le code implémente un algorithme EM pour un mélange de deux régresseurs (il devrait être relativement facile de l'étendre à des mélanges plus volumineux). Le code semble être assez robuste pour les ensembles de données aléatoires. Toutefois, contrairement à la régression linéaire, les modèles de mélange ont des objectifs non convexes. Par conséquent, pour un jeu de données réel, vous devrez peut-être exécuter quelques essais avec différents points de départ aléatoires.

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()
utilisateur1149913
la source
25

Ailleurs dans ce fil, user1149913 fournit d'excellents conseils (définir un modèle probabiliste) et du code pour une approche puissante (estimation EM). Deux problèmes restent à résoudre:

  1. Comment gérer les écarts par rapport au modèle probabiliste (ce qui est très évident dans les données de 2011-2012 et quelque peu évident dans les ondulations des points moins inclinés).

  2. Comment identifier de bonnes valeurs de départ pour l'algorithme EM (ou tout autre algorithme).

Pour résoudre le problème n ° 2, utilisez une transformation de Hough . Il s'agit d'un algorithme de détection de caractéristiques qui, pour rechercher des étendues linéaires de caractéristiques, peut être efficacement calculé comme une transformation de Radon .

XyX,ydans la transformation de Hough. Lorsque les entités du graphe d'origine suivent une ligne commune ou se rapprochent suffisamment d'une ligne, les ensembles de courbes qu'elles produisent dans la transformation de Hough tendent à avoir une intersection commune correspondant à cette ligne. En trouvant ces points de plus grande intensité dans la transformation de Hough, nous pouvons lire de bonnes solutions au problème initial.

Pour commencer à utiliser ces données, j'ai d'abord découpé les éléments auxiliaires (axes, graduations et étiquettes) et, pour faire bonne mesure, découpé les points manifestement périphériques en bas à droite et saupoudrés le long de l'axe du bas. (Lorsque cette opération n'est pas supprimée, la procédure fonctionne toujours correctement, mais elle détecte également les axes, les images, les séquences linéaires de graduations, les séquences linéaires d'étiquettes et même les points situés sporadiquement sur l'axe inférieur!)

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

(Ceci et le reste du code sont dans Mathematica .)

Image recadrée

A chaque point de cette image correspond une plage de courbes étroite dans la transformation de Hough, visible ici. Ce sont des ondes sinusoïdales:

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

Transformation de Hough

Cela montre visuellement que la question est un problème de regroupement de lignes : la transformation de Hough la réduit à un problème de regroupement de points auquel nous pouvons appliquer n'importe quelle méthode de regroupement.

Dans ce cas, le regroupement est tellement clair qu'un simple post-traitement de la transformation de Hough a suffi. Pour identifier les emplacements de plus grande intensité dans la transformation, j'ai augmenté le contraste et rendu floue la transformation sur un rayon d'environ 1%: cela est comparable aux diamètres des points de tracé de l'image d'origine.

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

Transformation floue

Seuil du résultat réduit à deux blobs minuscules dont les centroïdes identifient raisonnablement les points de plus grande intensité: ceux-ci estiment les lignes ajustées.

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

0,777

Transformation binarisée à seuil

Le côté gauche de l'image correspond à une direction de 0 degré (horizontal) et, en regardant de gauche à droite, cet angle augmente linéairement jusqu'à 180 degrés. En interpolant, je calcule que les deux blobs sont centrés à 19 et 57,1 degrés, respectivement. Nous pouvons également lire les interceptions à partir des positions verticales des blobs. Cette information donne les ajustements initiaux:

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 la même manière, on peut calculer les intersections correspondant à ces pentes, en donnant ces ajustements:

Lignes ajustées

(La ligne rouge correspond au minuscule point rose de la photo précédente et la ligne bleue à la plus grande goutte aqua.)

Dans une large mesure, cette approche a automatiquement résolu le premier problème: les écarts par rapport à la linéarité effacent les points d'intensité maximale, mais ne les modifient généralement pas beaucoup. Franchement, les points périphériques contribueront au bruit de bas niveau tout au long de la transformation de Hough, qui disparaîtra au cours des procédures de post-traitement.

À ce stade, vous pouvez fournir ces estimations en tant que valeurs de départ pour l'algorithme EM ou pour un minimiseur de vraisemblance (qui, avec de bonnes estimations, convergeront rapidement). Cependant, il serait préférable d’utiliser un estimateur de régression robuste tel que les moindres carrés repondérés de manière itérative . Il est capable de fournir un poids de régression à chaque point. Des poids faibles indiquent qu'un point n'appartient pas à une ligne. Exploitez ces poids, si vous le souhaitez, pour affecter chaque point à la ligne appropriée. Ensuite, après avoir classé les points, vous pouvez utiliser les moindres carrés ordinaires (ou toute autre procédure de régression) séparément sur les deux groupes de points.

whuber
la source
1
Les images disent mille mots et vous en avez 5. C’est un travail incroyable à partir d’un graphique rapide que j’ai fait juste pour les besoins de cette question! Gloire!
Jbbiomed
2
La transformation de Hough est largement utilisée dans le domaine de la vision par ordinateur pour identifier des lignes droites dans une image. Pourquoi ne devrait-il pas être utilisé dans les statistiques également? ;)
Lucas Reis
Xy
Oui. Imaginez, par exemple, le nombre de valeurs aberrantes impliquées dans la comparaison de deux images afin de déterminer si elles proviennent du même sujet. Et surtout, imaginez avoir à le faire en temps réel. La "rapidité" est un facteur très important dans la vision par ordinateur, et pas si important dans la statistique.
Lucas Reis
@RoyalTS Merci d'avoir souligné la nécessité de corriger l'un des extraits de code. Au moment où j'ai trouvé votre modification suggérée, elle avait été rejetée (à juste titre, parce que ce n'était pas tout à fait correct, mais cela ne fait rien: je vous suis reconnaissant d'avoir constaté une erreur). Je l'ai corrigé en supprimant la référence à rotation, qui avait été fixée à zéro à l'origine et ne faisait donc aucune différence.
whuber
15

J'ai trouvé cette question liée à une autre question . J'ai en fait fait des recherches universitaires sur ce genre de problème. S'il vous plaît vérifier ma réponse "racine carrée" montage? Une méthode d'adaptation avec plusieurs minima pour plus de détails.

L'approche basée sur la transformation de Hough de whuber est une très bonne solution pour les scénarios simples comme celui que vous avez donné. J'ai travaillé sur des scénarios avec des données plus complexes, telles que:

problème d'association de données - ensemble de données de bonbons

Mes co-auteurs et moi avons noté ce problème comme un problème d'association de données. Lorsque vous essayez de le résoudre, le problème principal est généralement combinatoire en raison de la quantité exponentielle de combinaisons de données possibles.

Nous avons une publication " Mélanges superposés de processus gaussiens pour le problème d'association de données " où nous avons abordé le problème général des courbes à N avec une technique itérative, donnant de très bons résultats. Vous pouvez trouver le code Matlab lié dans le document.

[Mise à jour] Une implémentation Python de la technique OMGP est disponible dans la bibliothèque GPClust .

J'ai un autre article dans lequel nous avons assoupli le problème afin d'obtenir un problème d'optimisation convexe, mais celui-ci n'a pas encore été accepté pour publication. Il est spécifique à 2 courbes, donc cela fonctionnerait parfaitement sur vos données. Fais moi savoir si tu es intéressé.

Steven
la source
1
Je suis triste de constater que plus de deux ans plus personne n'a voté pour cette réponse originale et utile. Entre-temps, le dernier document que vous avez mentionné a-t-il été accepté?
whuber
1
Le papier a bien été accepté il y a quelques mois à peine. Vous pouvez le télécharger ici gtas.unican.es/pub/378 . C'est en fait un problème assez rare (ce qui peut expliquer son manque de popularité), mais nous avons tout de même réussi à trouver des applications intéressantes. Jetez un coup d'œil aux expériences à la fin du document si vous le souhaitez.
Steven
2

user1149913 a une excellente réponse (+1), mais il me semble que votre collection de données s'est effondrée à la fin de 2011; il vous faudrait donc supprimer cette partie de vos données, puis exécuter plusieurs fois des choses avec différents aléatoires. coefficients de départ pour voir ce que vous obtenez.

Une façon simple de faire les choses serait de séparer vos données en deux ensembles à l’œil, puis d’utiliser la technique de modèle linéaire à laquelle vous êtes habitué. Dans R, ce serait la lmfonction.

Ou adapter deux lignes à l'œil. En R, vous utiliseriez ablinecela.

Les données sont mélangées, ont des valeurs aberrantes et se défont à la fin. Pourtant, vous avez deux lignes assez évidentes. Par conséquent, je ne suis pas sûr qu'une méthode sophistiquée en vaille la peine.

Wayne
la source