Рассчитать корреляцию в xarray с отсутствующими данными

Я пытаюсь рассчитать корреляцию между двумя наборами данных в xarray по временному измерению. Мой набор данных - это широта x долгота x время. В одном из моих наборов данных отсутствует достаточно данных, которые нецелесообразно интерполировать и устранять пробелы, вместо этого я хотел бы просто игнорировать отсутствующие значения. У меня есть несколько простых фрагментов кода, которые отчасти работают, но ни один из них не подходит для моего точного варианта использования. Например:

def covariance(x,y,dims=None):
    return xr.dot(x-x.mean(dims), y-y.mean(dims), dims=dims) / x.count(dims)

def correlation(x,y,dims=None):
    return covariance(x,y,dims) / (x.std(dims) * y.std(dims))

работает хорошо, если данные отсутствуют, но, конечно, не может работать с nans. Хотя есть хороший пример , написанный здесь для xarray, даже с этим кодом я изо всех сил пытаюсь вычислить корреляцию Пирсона, а не копейщика.

import numpy as np
import xarray as xr
import bottleneck

def covariance_gufunc(x, y):
    return ((x - x.mean(axis=-1, keepdims=True))
            * (y - y.mean(axis=-1, keepdims=True))).mean(axis=-1)

def pearson_correlation_gufunc(x, y):
    return covariance_gufunc(x, y) / (x.std(axis=-1) * y.std(axis=-1))

def spearman_correlation_gufunc(x, y):
    x_ranks = bottleneck.rankdata(x, axis=-1)
    y_ranks = bottleneck.rankdata(y, axis=-1)
    return pearson_correlation_gufunc(x_ranks, y_ranks)

def spearman_correlation(x, y, dim):
    return xr.apply_ufunc(
        spearman_correlation_gufunc, x, y,
        input_core_dims=[[dim], [dim]],
        dask='parallelized',
        output_dtypes=[float])

Наконец, было полезное обсуждение на github добавления этой функции в xarray, но она еще предстоит реализовать. Есть ли эффективный способ сделать это для наборов данных с пробелами в данных?


python numpy correlation python-xarray
person clifgray    schedule 24.10.2019    source источник

Ответы (2)


arrow_upward
4
arrow_downward

Я слежу за этим обсуждением Github и последующими попытками реализовать метод .corr(), кажется, мы довольно близки, но его еще нет.

Между тем, основной код, который большинство пытается объединить, довольно хорошо описан в этом другом ответе (Как применить линейную регрессию к каждому пикселю в большом многомерном массиве, содержащем NaN?). Это хорошее решение, которое использует векторизованные операции в NumPy и с небольшой настройкой (см. Принятый ответ в ссылке) может быть сделано для учета NaN по оси времени.

def lag_linregress_3D(x, y, lagx=0, lagy=0):
"""
Input: Two xr.Datarrays of any dimensions with the first dim being time. 
Thus the input data could be a 1D time series, or for example, have three 
dimensions (time,lat,lon). 
Datasets can be provided in any order, but note that the regression slope 
and intercept will be calculated for y with respect to x.
Output: Covariance, correlation, regression slope and intercept, p-value, 
and standard error on regression between the two datasets along their 
aligned time dimension.  
Lag values can be assigned to either of the data, with lagx shifting x, and
lagy shifting y, with the specified lag amount. 
""" 
#1. Ensure that the data are properly alinged to each other. 
x,y = xr.align(x,y)

#2. Add lag information if any, and shift the data accordingly
if lagx!=0:

    # If x lags y by 1, x must be shifted 1 step backwards. 
    # But as the 'zero-th' value is nonexistant, xr assigns it as invalid 
    # (nan). Hence it needs to be dropped
    x   = x.shift(time = -lagx).dropna(dim='time')

    # Next important step is to re-align the two datasets so that y adjusts
    # to the changed coordinates of x
    x,y = xr.align(x,y)

if lagy!=0:
    y   = y.shift(time = -lagy).dropna(dim='time')
    x,y = xr.align(x,y)

#3. Compute data length, mean and standard deviation along time axis: 
n = y.notnull().sum(dim='time')
xmean = x.mean(axis=0)
ymean = y.mean(axis=0)
xstd  = x.std(axis=0)
ystd  = y.std(axis=0)

#4. Compute covariance along time axis
cov   =  np.sum((x - xmean)*(y - ymean), axis=0)/(n)

#5. Compute correlation along time axis
cor   = cov/(xstd*ystd)

#6. Compute regression slope and intercept:
slope     = cov/(xstd**2)
intercept = ymean - xmean*slope  

#7. Compute P-value and standard error
#Compute t-statistics
tstats = cor*np.sqrt(n-2)/np.sqrt(1-cor**2)
stderr = slope/tstats

from scipy.stats import t
pval   = t.sf(tstats, n-2)*2
pval   = xr.DataArray(pval, dims=cor.dims, coords=cor.coords)

return cov,cor,slope,intercept,pval,stderr

Надеюсь это поможет! Скрестим пальцы, скоро для этого наступит слияние.

person AWilliams3142    schedule 22.02.2020
comment
PR для xr.cov и xr.corr был объединен на github.com/pydata/xarray/pull/4089, должен появиться в грядущей версии xarray v0.16.0. - person weiji14; 28.05.2020
comment
xr.corr и xr.cov уже выпущены! :) - person AWilliams3142; 14.07.2020

arrow_upward
0
arrow_downward

Решение находится в потоке github https://github.com/pydata/xarray/issues/1115

def covariance(x, y, dim=None):
    valid_values = x.notnull() & y.notnull()
    valid_count = valid_values.sum(dim)

    demeaned_x = (x - x.mean(dim)).fillna(0)
    demeaned_y = (y - y.mean(dim)).fillna(0)
    
    return xr.dot(demeaned_x, demeaned_y, dims=dim) / valid_count

def correlation(x, y, dim=None):
    # dim should default to the intersection of x.dims and y.dims
    return covariance(x, y, dim) / (x.std(dim) * y.std(dim))
person Zakari    schedule 31.07.2020