Найти матрицу корреляции

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

Просто быстро проверяя результаты на случайных данных, он возвращает то же самое, что и приведенный выше код @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((x-y)).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 должен делать именно то, что вы хотите, и это намного быстрее.

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

Посмотрев на ответ Джо Кингтона, я решил изучить код 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 () , потому что при этом не вычисляются избыточные коэффициенты и некоторые другие вещи. Это также будет быстрее и не вызовет проблем с памятью с матрицей из 50 000 строк, потому что тогда вы можете либо сохранить каждый набор r, либо обработать их перед созданием другого набора. Не сохраняя никаких долгосрочных значений r, я смог заставить приведенный выше код работать на наборе случайным образом сгенерированных данных размером 50 000 x 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-значением - 8,8x, если я оставлю там материал с p-значением (я использовал 10 столбцов с сотнями строк). Я также проверил, дает ли он те же результаты. Это небольшое улучшение, но оно может помочь.

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

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

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

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

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

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