python - corrcoef - Calcular el coeficiente de correlación entre dos matrices multidimensionales

corrcoef python

Question

Tengo dos matrices que tienen las formas NXT y MXT . Me gustaría calcular el coeficiente de correlación a través de T entre cada posible par de filas n y m (de N y M , respectivamente).

¿Cuál es la forma más rápida y más pitónica de hacer esto? (Dar vueltas sobre N y M me parece que no es ni rápido ni pitónico.) Espero que la respuesta implique numpy y / o scipy . En este momento, mis matrices son matrices numpy , pero estoy dispuesto a convertirlas a un tipo diferente.

Espero que mi salida sea una matriz con la forma NXM .

Nota: cuando digo "coeficiente de correlación", me refiero al coeficiente de correlación momento-producto de Pearson .

Aquí hay algunas cosas a tener en cuenta:

• La correlate función numpy requiere que las matrices de entrada sean unidimensionales.
• La función corrcoef acepta matrices bidimensionales, pero deben tener la misma forma.
• La función pearsonr requiere que las matrices de entrada sean unidimensionales.

Answer 1

@Divakar ofrece una gran opción para calcular la correlación sin escala, que es lo que pedí originalmente.

Para calcular el coeficiente de correlación, se requiere un poco más:

import numpy as np

def generate_correlation_map(x, y):
    """Correlate each n with each m.

    Parameters
    ----------
    x : np.array
      Shape N X T.

    y : np.array
      Shape M X T.

    Returns
    -------
    np.array
      N X M array in which each element is a correlation coefficient.

    """
    mu_x = x.mean(1)
    mu_y = y.mean(1)
    n = x.shape[1]
    if n != y.shape[1]:
        raise ValueError(''x and y must '' +
                         ''have the same number of timepoints.'')
    s_x = x.std(1, ddof=n - 1)
    s_y = y.std(1, ddof=n - 1)
    cov = np.dot(x,
                 y.T) - n * np.dot(mu_x[:, np.newaxis],
                                  mu_y[np.newaxis, :])
    return cov / np.dot(s_x[:, np.newaxis], s_y[np.newaxis, :])

Aquí hay una prueba de esta función, que pasa:

from scipy.stats import pearsonr

def test_generate_correlation_map():
    x = np.random.rand(10, 10)
    y = np.random.rand(20, 10)
    desired = np.empty((10, 20))
    for n in range(x.shape[0]):
        for m in range(y.shape[0]):
            desired[n, m] = pearsonr(x[n, :], y[m, :])[0]
    actual = generate_correlation_map(x, y)
    np.testing.assert_array_almost_equal(actual, desired)

Answer 2

Correlación (caso "válido" predeterminado) entre dos matrices 2D:

Simplemente puede usar la matriz de multiplicación np.dot así -

out = np.dot(arr_one,arr_two.T)

La correlación con el caso "valid" predeterminado entre cada combinación de filas por pares (fila1, fila2) de las dos matrices de entrada correspondería al resultado de multiplicación en cada posición (fila1, fila2).

Cálculo del coeficiente de correlación en filas para dos matrices 2D:

def corr2_coeff(A,B):
    # Rowwise mean of input arrays & subtract from input arrays themeselves
    A_mA = A - A.mean(1)[:,None]
    B_mB = B - B.mean(1)[:,None]

    # Sum of squares across rows
    ssA = (A_mA**2).sum(1);
    ssB = (B_mB**2).sum(1);

    # Finally get corr coeff
    return np.dot(A_mA,B_mB.T)/np.sqrt(np.dot(ssA[:,None],ssB[None]))

Esto se basa en esta solución de How to apply corr2 functions in Multidimentional arrays in MATLAB

Benchmarking

Esta sección compara el rendimiento en tiempo de ejecución con el enfoque propuesto con el enfoque basado en generate_correlation_map y pearsonr se enumera en la otra respuesta. (tomado de la función test_generate_correlation_map() sin el código de verificación de corrección de valores al final de la misma). Tenga en cuenta que los tiempos para el enfoque propuesto también incluyen una verificación al comienzo para verificar el mismo número de columnas en las dos matrices de entrada, como también se hizo en esa otra respuesta. Los tiempos de ejecución se enumeran a continuación.

Caso 1:

In [106]: A = np.random.rand(1000,100)

In [107]: B = np.random.rand(1000,100)

In [108]: %timeit corr2_coeff(A,B)
100 loops, best of 3: 15 ms per loop

In [109]: %timeit generate_correlation_map(A, B)
100 loops, best of 3: 19.6 ms per loop

Caso # 2:

In [110]: A = np.random.rand(5000,100)

In [111]: B = np.random.rand(5000,100)

In [112]: %timeit corr2_coeff(A,B)
1 loops, best of 3: 368 ms per loop

In [113]: %timeit generate_correlation_map(A, B)
1 loops, best of 3: 493 ms per loop

Caso # 3:

In [114]: A = np.random.rand(10000,10)

In [115]: B = np.random.rand(10000,10)

In [116]: %timeit corr2_coeff(A,B)
1 loops, best of 3: 1.29 s per loop

In [117]: %timeit generate_correlation_map(A, B)
1 loops, best of 3: 1.83 s per loop

El otro enfoque pearsonr based parecía demasiado lento, pero aquí están los tiempos de ejecución para un pequeño tamaño de datos:

In [118]: A = np.random.rand(1000,100)

In [119]: B = np.random.rand(1000,100)

In [120]: %timeit corr2_coeff(A,B)
100 loops, best of 3: 15.3 ms per loop

In [121]: %timeit generate_correlation_map(A, B)
100 loops, best of 3: 19.7 ms per loop

In [122]: %timeit pearsonr_based(A,B)
1 loops, best of 3: 33 s per loop