Поиск корреляционной матрицы

Новое решение

Посмотрев на ответ Джо Кингтона, я решил изучить код corrcoef() и был вдохновлен им сделать следующую реализацию.

 ms = data.mean(axis=1)[(slice(None,None,None),None)] datam = data - ms datass = np.sqrt(scipy.stats.ss(datam,axis=1)) for i in xrange(rows): temp = np.dot(datam[i:],datam[i].T) rs = temp / (datass[i:]*datass[i]) 

Каждый цикл через генерирует коэффициенты Пирсона между строками i и строками i до последней строки. Это очень быстро. Он по крайней мере в 1,5 corrcoef() быстрее, чем использование corrcoef() потому что он не избыточно вычисляет коэффициенты и несколько других вещей. Он также будет быстрее и не даст вам проблем с памятью с матрицей строк в 50 000, потому что тогда вы можете либо сохранить каждый набор r, либо обработать их до создания другого набора. Не сохраняя ни одного из долговременных сроков, я смог получить приведенный выше код для запуска на 50 000 х 10 наборов случайно генерируемых данных за минуту на моем довольно новом ноутбуке.

Старое решение

Во-первых, я бы не рекомендовал распечатывать r на экране. Для 100 строк (10 столбцов) это разница в 19,79 секунды с печатью против 0,301 секунды без использования вашего кода. Просто сохраните r и используйте их позже, если хотите, или выполните некоторую обработку на них, когда вы идете вперед, ища некоторые из самых больших r.

Во-вторых, вы можете получить некоторую экономию, не избыточно вычисляя некоторые количества. Коэффициент Пирсона рассчитывается в scipy с использованием некоторых величин, которые вы можете предварительно вычислять, а не вычислять каждый раз, когда используется строка. Кроме того, вы не используете p-значение (которое также возвращается pearsonr() поэтому давайте поцарапаем это тоже. Используя приведенный ниже код:

 r = np.zeros((rows,rows)) ms = data.mean(axis=1) datam = np.zeros_like(data) for i in xrange(rows): datam[i] = data[i] - ms[i] datass = scipy.stats.ss(datam,axis=1) for i in xrange(rows): for j in xrange(i,rows): r_num = np.add.reduce(datam[i]*datam[j]) r_den = np.sqrt(datass[i]*datass[j]) r[i,j] = min((r_num / r_den), 1.0) 

Я получаю ускорение около 4.8x по сравнению с прямым scipy-кодом, когда я удалял материал p-value – 8.8x, если я оставил там p-значение (я использовал 10 столбцов с сотнями строк). Я также проверил, что он дает те же результаты. Это не очень большое улучшение, но это может помочь.

В конечном счете, вы столкнулись с проблемой, которую вы вычисляете (50000) * (50001) / 2 = 1 250 025 000 коэффициентов Пирсона (если я правильно рассчитываю). Это много. Кстати, нет необходимости вычислять коэффициент Pearson каждой строки с самим собой (он будет равен 1), но это только избавит вас от вычисления 50 000 коэффициентов Пирсона. С приведенным выше кодом я ожидаю, что для вычисления будет потребоваться около 4 1/4 часа, если у вас есть 10 столбцов для ваших данных, основанных на моих результатах на меньших наборах данных.

Вы можете получить некоторое улучшение, взяв вышеуказанный код в Cython или что-то подобное. Я ожидаю, что вы, возможно, получите 10-кратное улучшение по сравнению с прямым Scipy, если вам повезет. Кроме того, как предложено pyInTheSky, вы можете выполнить некоторую многопроцессорную обработку.

Вы пробовали использовать numpy.corrcoef ? Видя, как вы не используете p-значения, он должен делать именно то, что вы хотите, с наименьшими проблемами. (Если я не помню точно, что такое pearson's R, что вполне возможно.)

Просто быстро проверяя результаты по случайным данным, он возвращает точно такую ​​же вещь, как и код @Justin Peel выше и работает ~ 100 раз быстрее.

Например, тестирование объектов с 1000 строк и 10 столбцов случайных данных …:

 import numpy as np import scipy as sp import scipy.stats def main(): data = np.random.random((1000, 10)) x = corrcoef_test(data) y = justin_peel_test(data) print 'Maximum difference between the two results:', np.abs((xy)).max() return data def corrcoef_test(data): """Just using numpy's built-in function""" return np.corrcoef(data) def justin_peel_test(data): """Justin Peel's suggestion above""" rows = data.shape[0] r = np.zeros((rows,rows)) ms = data.mean(axis=1) datam = np.zeros_like(data) for i in xrange(rows): datam[i] = data[i] - ms[i] datass = sp.stats.ss(datam,axis=1) for i in xrange(rows): for j in xrange(i,rows): r_num = np.add.reduce(datam[i]*datam[j]) r_den = np.sqrt(datass[i]*datass[j]) r[i,j] = min((r_num / r_den), 1.0) r[j,i] = r[i,j] return r data = main() 

Достигает максимальной абсолютной разницы ~ 3.3e-16 между двумя результатами

И тайминги:

 In [44]: %timeit corrcoef_test(data) 10 loops, best of 3: 71.7 ms per loop In [45]: %timeit justin_peel_test(data) 1 loops, best of 3: 6.5 s per loop 

numpy.corrcoef должен делать то, что вы хотите, и это намного быстрее.

вы можете использовать модуль мультипроцесса python, разбивать строки на 10 наборов, буферизировать результаты и затем печатать материал (это только ускорит его на многоядерной машине)

http://docs.python.org/library/multiprocessing.html

Кстати, вам также придется превратить ваш фрагмент в функцию, а также подумать о том, как выполнить сборку данных. у каждого подпроцесса есть список, подобный этому … [startcord, stopcord, buff] .. может работать красиво

 def myfunc(thelist): for i in xrange(thelist[0]:thelist[1]): .... thelist[2] = result