Os dados têm duas tendências; como extrair linhas de tendência independentes?

34

Eu tenho um conjunto de dados que não são ordenados de maneira específica, mas quando plotados claramente têm duas tendências distintas. Uma regressão linear simples não seria realmente adequada aqui devido à clara distinção entre as duas séries. Existe uma maneira simples de obter as duas linhas de tendência lineares independentes?

Para constar, estou usando Python e estou razoavelmente confortável com programação e análise de dados, incluindo aprendizado de máquina, mas estou disposto a pular para R se for absolutamente necessário.

insira a descrição da imagem aqui

jbbiomed
fonte
6
Melhor resposta que eu tenho até agora é para imprimir esta no papel de gráfico e use um lápis e régua e calculadora ...
jbbiomed
Talvez você possa calcular declives entre pares e agrupá-los em dois "agrupamentos de declive". No entanto, isso falhará se você tiver duas tendências paralelas.
Thomas Jungblut
11
Não tenho nenhuma experiência pessoal com isso, mas acho que vale a pena conferir os modelos de estatísticas . Estatisticamente, uma regressão linear com uma interação de grupo seria adequada (a menos que você está dizendo que você tem dados não agrupadas, caso em que isso é um pouco mais peludo ...)
Matt Parker
11
Infelizmente, não se trata de dados de efeito, mas de uso, e claramente do uso de dois sistemas separados, misturados no mesmo conjunto de dados. Quero ser capaz de descrever os dois padrões de uso, mas não posso voltar atrás e relembrar dados, pois isso representa cerca de 6 anos de informações coletadas por um cliente.
Jbbiomed
2
Apenas para garantir: seu cliente não possui dados adicionais que indiquem quais medidas são provenientes de qual população? Isso é 100% dos dados que você ou seu cliente tem ou pode encontrar. Além disso, 2012 parece que sua coleta de dados desmoronou ou um ou ambos os sistemas caíram no chão. Me faz pensar se as tendências até esse ponto importam muito.
Wayne

Respostas:

30

Para resolver seu problema, uma boa abordagem é definir um modelo probabilístico que corresponda às suposições sobre o seu conjunto de dados. No seu caso, você provavelmente deseja uma mistura de modelos de regressão linear. Você pode criar um modelo de "mistura de regressores" semelhante a um modelo de mistura gaussiano associando diferentes pontos de dados a diferentes componentes de mistura.

Eu incluí algum código para você começar. O código implementa um algoritmo EM para uma mistura de dois regressores (deve ser relativamente fácil estender para misturas maiores). O código parece ser bastante robusto para conjuntos de dados aleatórios. No entanto, diferentemente da regressão linear, os modelos de mistura têm objetivos não convexos; portanto, para um conjunto de dados real, você pode precisar executar algumas tentativas com diferentes pontos de partida aleatórios.

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()
user1149913
fonte
25

Em outras partes deste segmento, o usuário1149913 fornece ótimos conselhos (define um modelo probabilístico) e codifica uma abordagem poderosa (estimativa EM). Ainda há duas questões a serem abordadas:

  1. Como lidar com desvios do modelo probabilístico (que são muito evidentes nos dados de 2011-2012 e um tanto evidentes nas ondulações dos pontos menos inclinados).

  2. Como identificar bons valores iniciais para o algoritmo EM (ou qualquer outro algoritmo).

Para resolver o problema nº 2, considere usar uma transformação Hough . Este é um algoritmo de detecção de recurso que, para encontrar trechos lineares de recursos, pode ser computado com eficiência como uma transformação de Radon .

xyx,yna transformação Hough. Quando os recursos no gráfico original caem ao longo de uma linha comum, ou quase o suficiente para um, as coleções de curvas que eles produzem na transformação Hough tendem a ter uma interseção comum correspondente a essa linha comum. Ao encontrar esses pontos de maior intensidade na transformação de Hough, podemos ler boas soluções para o problema original.

Para começar com esses dados, recortei primeiro o material auxiliar (eixos, marcas de escala e rótulos) e, na medida certa, recortei os pontos obviamente periféricos no canto inferior direito e salpiquei ao longo do eixo inferior. (Quando esse material não é cortado, o procedimento ainda funciona bem, mas também detecta os eixos, os quadros, as seqüências lineares de ticks, as sequências lineares de etiquetas e até os pontos esporadicamente no eixo 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 e o restante do código estão no Mathematica .)

Imagem recortada

Para cada ponto nesta imagem corresponde uma faixa estreita de curvas na transformação Hough, visível aqui. São ondas senoidais:

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

Transformação Hough

Isso torna visualmente manifesto o sentido em que a questão é um problema de agrupamento de linhas : a transformação Hough a reduz a um problema de agrupamento de pontos , ao qual podemos aplicar qualquer método de agrupamento que desejar.

Nesse caso, o armazenamento em cluster é tão claro que o pós-processamento simples da transformação Hough foi suficiente. Para identificar os locais de maior intensidade na transformação, aumentei o contraste e embaçei a transformação em um raio de cerca de 1%: é comparável aos diâmetros dos pontos de plotagem na imagem original.

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

Transformação turva

Limiar o resultado reduziu-o a duas pequenas bolhas cujos centróides identificam razoavelmente os pontos de maior intensidade: estimam as linhas ajustadas.

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

0,777

Transformação binarizada com limite

O lado esquerdo da imagem corresponde a uma direção de 0 graus (horizontal) e, conforme olhamos da esquerda para a direita, esse ângulo aumenta linearmente para 180 graus. Interpolando, calculo que os dois blobs estão centrados em 19 e 57,1 graus, respectivamente. Também podemos ler as interceptações nas posições verticais dos blobs. Esta informação produz os ajustes iniciais:

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 maneira semelhante, é possível calcular as interceptações correspondentes a essas pistas, dando os seguintes ajustes:

Linhas ajustadas

(A linha vermelha corresponde ao pequeno ponto rosa na imagem anterior e a linha azul corresponde ao blob aqua maior.)

Em grande parte, essa abordagem lidou automaticamente com a primeira questão: os desvios da linearidade eliminam os pontos de maior intensidade, mas normalmente não os alteram muito. Pontos francamente afastados contribuirão com ruído de baixo nível em toda a transformação Hough, que desaparecerá durante os procedimentos de pós-processamento.

Nesse ponto, é possível fornecer essas estimativas como valores iniciais para o algoritmo EM ou como um minimizador de probabilidade (que, dadas boas estimativas, convergirá rapidamente). Melhor, porém, seria usar um estimador de regressão robusto, como mínimos quadrados ponderados iterativamente . É capaz de fornecer um peso de regressão para cada ponto. Pesos baixos indicam que um ponto não "pertence" a uma linha. Explore esses pesos, se desejado, para atribuir cada ponto à sua linha correta. Depois de classificar os pontos, você pode usar mínimos quadrados comuns (ou qualquer outro procedimento de regressão) separadamente nos dois grupos de pontos.

whuber
fonte
11
As imagens dizem mais que mil palavras e você tem 5. Este é um trabalho incrível de um gráfico rápido que fiz apenas para o propósito desta pergunta! Parabéns!
Jbbiomed
2
A transformação de Hough é amplamente usada no campo Computer Vision para identificar linhas retas em uma imagem. Por que não deveria ser usado também nas estatísticas? ;)
Lucas Reis
xy
Sim. Imagine, por exemplo, a quantidade de discrepantes envolvidos na comparação de duas imagens para detectar se são do mesmo assunto. E, acima de tudo, imagine ter que fazer isso em tempo real. "Velocidade" é um fator muito importante na visão computacional, e não tão importante na estatística.
Lucas Reis
@RoyalTS Obrigado por apontar a necessidade de uma correção para um dos trechos de código. No momento em que encontrei sua alteração sugerida, ela havia sido rejeitada (corretamente, porque não estava certa, mas não se esqueça disso: agradeço que você tenha notado um erro). Corrigi-o removendo a referência a rotation, que foi definida como zero originalmente e, portanto, não fez diferença.
whuber
15

Encontrei esta questão vinculada a outra questão . Na verdade, eu fiz pesquisas acadêmicas sobre esse tipo de problema. Por favor, verifique a minha resposta "Menos raiz quadrada". Um método de ajuste com vários mínimos para obter mais detalhes.

A abordagem baseada em transformação de Hough da whuber's é uma solução muito boa para cenários simples como o que você deu. Trabalhei em cenários com dados mais complexos, como este:

problema de associação de dados - conjunto de dados candy

Meus co-autores e eu denotamos isso como um problema de "associação de dados". Quando você tenta resolvê-lo, o principal problema geralmente é combinatório devido à quantidade exponencial de combinações de dados possíveis.

Temos uma publicação " Misturas sobrepostas de processos gaussianos para o problema de associação de dados ", onde abordamos o problema geral das curvas N com uma técnica iterativa, obtendo resultados muito bons. Você pode encontrar o código Matlab vinculado no documento.

[Atualização] Uma implementação em Python da técnica OMGP pode ser encontrada na biblioteca GPClust .

Tenho outro artigo em que relaxamos o problema para obter um problema de otimização convexa, mas ele ainda não foi aceito para publicação. É específico para 2 curvas, por isso funcionaria perfeitamente em seus dados. Deixe-me saber se você está interessado.

Steven
fonte
11
Fico triste ao ver que, ao longo de dois anos, ninguém mais votou nessa resposta original e valiosa. Enquanto isso, o último artigo que você mencionou foi aceito?
whuber
11
O documento foi de fato aceito, apenas alguns meses atrás. Você pode baixá-lo aqui gtas.unican.es/pub/378 . Na verdade, esse é um problema bastante raro (que pode explicar sua falta de popularidade), mas ainda conseguimos encontrar alguns aplicativos interessantes. Dê uma olhada nas experiências no final do artigo, se quiser.
1716 Steven
2

user1149913 tem uma excelente resposta (+1), mas parece-me que sua coleta de dados desmoronou no final de 2011; portanto, você teria que cortar essa parte dos dados e ainda executar algumas vezes com aleatoriamente diferentes coeficientes iniciais para ver o que você recebe.

Uma maneira simples de fazer as coisas seria separar seus dados em dois conjuntos por olho e, em seguida, usar qualquer técnica de modelo linear com a qual você esteja acostumado. Em R, seria a lmfunção.

Ou ajuste duas linhas a olho nu. Em R você usaria ablinepara fazer isso.

Os dados são confusos, apresentam discrepâncias e desmoronam no final, mas os olhos têm duas linhas bastante óbvias, então não tenho certeza se um método sofisticado vale a pena.

Wayne
fonte