Данные имеют две тенденции; как извлечь независимые линии тренда?

34

У меня есть набор данных, который не упорядочен каким-либо конкретным способом, но при четком графике имеет две четкие тенденции. Простая линейная регрессия здесь не совсем подходит из-за четкого различия между двумя рядами. Есть ли простой способ получить две независимые линейные линии тренда?

Для справки: я использую Python, и я достаточно хорошо разбираюсь в программировании и анализе данных, включая машинное обучение, но готов перейти на R, если это абсолютно необходимо.

введите описание изображения здесь

jbbiomed
источник
6
Лучший ответ, который у меня есть, это напечатать это на миллиметровой бумаге и использовать карандаш, линейку и калькулятор ...
jbbiomed
Может быть, вы можете вычислить попарные уклоны и сгруппировать их в два «склона-кластера». Однако это не удастся, если у вас есть две параллельные тенденции.
Томас Юнгблут
1
У меня нет никакого личного опыта с этим, но я думаю, что statsmodels стоило бы проверить. Статистически, линейная регрессия с взаимодействием для группы была бы адекватной (если только вы не говорите, что у вас есть разгруппированные данные, в этом случае это немного более опасно ...)
Мэтт Паркер
1
К сожалению, это не данные об эффектах, а данные об использовании, и, очевидно, данные об использовании из двух разных систем смешаны в одном наборе данных. Я хочу иметь возможность описать две модели использования, но я не могу вернуться и вспомнить данные, поскольку они представляют собой информацию, собранную клиентом за 6 лет.
jbbiomed
2
Просто чтобы убедиться: у вашего клиента нет никаких дополнительных данных, которые бы указывали, какие измерения получены от какой совокупности? Это 100% данных, которые вы или ваш клиент имеете или можете найти. Кроме того, 2012 год выглядит так, будто ваш сбор данных развалился или одна или обе ваши системы провалились. Заставляет меня задуматься, имеет ли значение линия тренда до этого момента.
Уэйн

Ответы:

30

Чтобы решить вашу проблему, хорошим подходом является определение вероятностной модели, которая соответствует предположениям о вашем наборе данных. В вашем случае вам, вероятно, понадобится смесь моделей линейной регрессии. Вы можете создать модель «смеси регрессоров», аналогичную модели смеси Гаусса, связав разные точки данных с различными компонентами смеси.

Я включил некоторый код, чтобы вы начали. Код реализует EM-алгоритм для смеси двух регрессоров (его должно быть относительно легко распространить на большие смеси). Код кажется достаточно надежным для случайных наборов данных. Однако, в отличие от линейной регрессии, смешанные модели имеют невыпуклые цели, поэтому для реального набора данных может потребоваться провести несколько испытаний с разными случайными начальными точками.

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
источник
25

В другом месте в этой теме user1149913 предоставляет полезные советы (определяют вероятностную модель) и код для мощного подхода (оценка EM). Осталось решить две проблемы:

  1. Как справиться с отклонениями от вероятностной модели (что очень очевидно в данных за 2011-2012 гг. И несколько очевидно в волнистости менее наклонных точек).

  2. Как определить хорошие начальные значения для алгоритма EM (или любого другого алгоритма).

Чтобы обратиться к № 2, подумайте об использовании преобразования Хафа . Это алгоритм обнаружения признаков, который для нахождения линейных участков признаков может быть эффективно вычислен как преобразование Радона .

xyx,yв преобразовании Хафа. Когда объекты в исходном графике располагаются вдоль общей линии или достаточно близко к одной, тогда наборы кривых, которые они создают в преобразовании Хафа, имеют общее пересечение, соответствующее этой общей линии. Найдя эти точки наибольшей интенсивности в преобразовании Хафа, мы можем зачитать хорошие решения исходной задачи.

Чтобы начать работу с этими данными, я сначала обрезал вспомогательные элементы (оси, отметки и метки) и для хорошей цели обрезал явно отдаленные точки внизу справа и посыпал вдоль нижней оси. (Когда этот материал не обрезается, процедура все еще работает хорошо, но она также обнаруживает оси, кадры, линейные последовательности тиков, линейные последовательности меток и даже точки, время от времени лежащие на нижней оси!)

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

(Этот и остальной код приведены в Mathematica .)

Обрезанное изображение

Каждой точке на этом изображении соответствует узкий диапазон кривых в преобразовании Хафа, видимый здесь. Это синусоиды:

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

Хау трансформация

Это наглядно демонстрирует смысл, в котором вопрос представляет собой проблему линейной кластеризации : преобразование Хафа сводит его к проблеме точечной кластеризации , к которой мы можем применить любой метод кластеризации, который нам нравится.

В этом случае кластеризация настолько ясна, что достаточно простой постобработки преобразования Хафа. Чтобы определить места наибольшей интенсивности в преобразовании, я увеличил контраст и размыл преобразование по радиусу около 1%: это сопоставимо с диаметрами точек графика на исходном изображении.

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

Затуманенное преобразование

Установление порога результата сузило его до двух крошечных сгустков, центроиды которых разумно идентифицируют точки наибольшей интенсивности: они оценивают подобранные линии.

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

0.777

Пороговое преобразование в двоичную форму

Левая сторона изображения соответствует направлению 0 градусов (по горизонтали), и, как мы видим слева направо, этот угол увеличивается линейно до 180 градусов. Интерполируя, я вычисляю, что эти две капли центрированы в 19 и 57,1 градусах соответственно. Мы также можем считывать перехваты с вертикальных позиций капель. Эта информация дает начальные соответствия:

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}

Аналогичным образом можно вычислить точки пересечения, соответствующие этим наклонам, давая следующие соответствия:

Установленные линии

(Красная линия соответствует крошечной розовой точке на предыдущем изображении, а синяя линия соответствует большему водному шарику.)

В значительной степени этот подход автоматически решает первую проблему: отклонения от линейности размывают точки наибольшей интенсивности, но обычно не сильно их смещают. Откровенно говоря, отдаленные точки будут вносить низкий уровень шума во всем преобразовании Хафа, который исчезнет во время процедур последующей обработки.

В этот момент можно предоставить эти оценки в качестве начальных значений для алгоритма EM или для минимизатора вероятности (который, при наличии хороших оценок, будет быстро сходиться). Однако лучше было бы использовать надежную регрессионную оценку, такую ​​как метод наименьших квадратов с повторным взвешиванием . Он способен обеспечить регрессионный вес для каждой точки. Низкие веса указывают, что точка не принадлежит линии. Используйте эти веса, если хотите, чтобы назначить каждую точку соответствующей линии. Затем, классифицировав точки, вы можете использовать обычные наименьшие квадраты (или любую другую процедуру регрессии) отдельно для двух групп точек.

Whuber
источник
1
Картинки рассказывают тысячу слов, а у вас их 5. Это невероятная работа из быстрого графика, который я сделал только для целей этого вопроса! Престижность!
jbbiomed
2
Преобразование Хафа широко используется в области компьютерного зрения для определения прямых линий на изображении. Почему его не следует использовать и в статистике? ;)
Лукас Рейс
xy
Да. Представьте, например, количество выбросов, участвующих в сравнении двух изображений, чтобы определить, относятся ли они к одному объекту. И, самое главное, представьте, что вам нужно делать это в реальном времени. «Скорость» является очень важным фактором в компьютерном зрении, и не столь важным в статистике.
Лукас Рейс
@RoyalTS Спасибо, что указали на необходимость исправления одного из фрагментов кода. К тому времени, как я нашел предложенное вами изменение, оно было отклонено (правильно, потому что оно было не совсем правильным, но не берите в голову: я рад, что вы заметили, что произошла ошибка). Я исправил это, удалив ссылку на rotation, которая изначально была установлена ​​на ноль и поэтому не имела никакого значения.
whuber
15

Я обнаружил, что этот вопрос связан с другим вопросом . Я действительно проводил академические исследования по этой проблеме. Пожалуйста, проверьте мой ответ "Наименьший квадратный корень" подходит? Подходящий метод с несколькими минимумами для более подробной информации.

Подход, основанный на преобразовании Хафа, является очень хорошим решением для простых сценариев, которые вы предложили. Я работал над сценариями с более сложными данными, такими как это:

проблема ассоциации данных - набор данных конфет

Мои соавторы и я обозначили это как проблему «ассоциации данных». Когда вы пытаетесь ее решить, основная проблема обычно является комбинаторной из-за экспоненциального количества возможных комбинаций данных.

У нас есть публикация « Перекрывающиеся смеси гауссовских процессов для задачи об ассоциации данных », в которой мы подошли к общей проблеме N кривых с помощью итерационной техники, дающей очень хорошие результаты. Вы можете найти код Matlab, связанный в статье.

[Обновление] Реализацию Python реализации технологии OMGP можно найти в библиотеке GPClust .

У меня есть еще одна статья, где мы ослабили задачу, чтобы получить задачу выпуклой оптимизации, но она еще не принята к публикации. Он специфичен для 2 кривых, поэтому он отлично подойдет для ваших данных. Дайте знать, если вас это заинтересовало.

Стивен
источник
1
Мне грустно видеть, что за два года никто больше не проголосовал за этот оригинальный и ценный ответ. Между тем, была ли принята последняя упомянутая вами статья?
whuber
1
Документ действительно был принят, всего несколько месяцев назад. Вы можете скачать его здесь gtas.unican.es/pub/378 . На самом деле это довольно редкая проблема (которая может объяснить его отсутствие популярности), но нам все же удалось найти несколько интересных приложений. Посмотрите на эксперименты в конце статьи, если хотите.
Стивен
2

У user1149913 отличный ответ (+1), но мне кажется, что ваш сбор данных развалился в конце 2011 года, так что вам пришлось бы отрезать эту часть своих данных, а затем несколько раз запускать вещи с разными случайными числами. Начальные коэффициенты, чтобы увидеть, что вы получите.

Одним из простых способов сделать это было бы разделить ваши данные на два набора на глаз, а затем использовать любой метод линейной модели, к которому вы привыкли. В R это была бы lmфункция.

Или подходят две линии на глаз. В R вы бы использовали ablineдля этого.

Данные перемешаны, имеют выбросы и в конце распадаются, но на глаз есть две довольно очевидные линии, поэтому я не уверен, что причудливый метод того стоит.

Wayne
источник