Как применить линейную регрессию к каждому пикселю в большом многомерном массиве, содержащем NaN?

Question

У меня есть одномерный массив значений независимых переменных (x_array), которые соответствуют временным шагам в трехмерном массиве пространственных данных с несколькими временными шагами (y_array). Мои фактические данные намного больше: 300+ временных шагов и до 3000 * 3000 пикселей:

import numpy as np
from scipy.stats import linregress

# Independent variable: four time-steps of 1-dimensional data 
x_array = np.array([0.5, 0.2, 0.4, 0.4])

# Dependent variable: four time-steps of 3x3 spatial data
y_array = np.array([[[-0.2,   -0.2,   -0.3],
                     [-0.3,   -0.2,   -0.3],
                     [-0.3,   -0.4,   -0.4]],

                    [[-0.2,   -0.2,   -0.4],
                     [-0.3,   np.nan, -0.3],
                     [-0.3,   -0.3,   -0.4]],

                    [[np.nan, np.nan, -0.3],
                     [-0.2,   -0.3,   -0.7],
                     [-0.3,   -0.3,   -0.3]],

                    [[-0.1,   -0.3,   np.nan],
                     [-0.2,   -0.3,   np.nan],
                     [-0.1,   np.nan, np.nan]]])

Я хочу вычислить линейную регрессию на пиксель и получить R-квадрат, P-значения, перехваты и наклоны для каждого xy пикселя в y_array, со значениями для каждого временного шага в x_array в качестве моей независимой переменной ,

Я могу изменить форму, чтобы получить данные в форме для ввода их в np.polyfit, который векторизован и быстр:

# Reshape so rows = number of time-steps and columns = pixels:
y_array_reshaped = y_array.reshape(len(y_array), -1)

# Do a first-degree polyfit
np.polyfit(x_array, y_array_reshaped, 1)

Однако это игнорирует пиксели, которые содержат любые значения NaN (np.polyfit не поддерживает значения NaN), и не вычисляет требуемую статистику (R-квадрат, P-значения, перехваты и уклоны).

Ответ здесь использует scipy.stats import linregress, который вычисляет необходимую мне статистику и предлагает избегать проблем NaN, маскируя эти NaN значения. Однако этот пример относится к двум одномерным массивам, и я не могу понять, как применить подобный подход маскирования к моему случаю, когда каждый столбец в y_array_reshaped будет иметь различный набор значений NaN.

Мой вопрос: Как рассчитать статистику регрессии для каждого пикселя в большом многомерном массиве (300 x 3000 x 3000), содержащем множество значений NaN, достаточно быстрым, векторизованным способом?

Желаемый результат: Массив статистических значений регрессии 3 x 3 (например, R-квадрат) для каждого пикселя в y_array, даже если этот пиксель содержит значения NaN в некоторый момент времени серия

Answer 1

Ответ, предоставленный здесь https://hrishichandanpurkar.blogspot.com/2017/09/vectorized-functions-for-correlation.html, абсолютно хорош в том смысле, что он в основном использует огромные возможности numpy векторизации и вещания, но предполагает, что анализируемые данные полны, что обычно не имеет место в реальных условиях. исследовательский цикл. Один ответ выше предназначен для решения проблемы отсутствующих данных, но я лично думаю, что нужно обновить больше кодов просто потому, что np.mean() вернет nan, если в данных присутствует nan. К счастью, numpy предоставил nanmean(), nanstd() и т. Д., Которые мы можем использовать для вычисления среднего значения, стандартной ошибки и т. Д., Игнорируя nans в данных. Между тем, программа в оригинальном блоге ориентирована на данные в формате netCDF. Некоторые могут этого не знать, но лучше знакомы с форматом numpy.array. Поэтому я приведу здесь пример кода, показывающий, как рассчитать ковариацию, коэффициенты корреляции и т. Д. Между двумя трехмерными массивами (размерность n-D имеет одинаковую логику). Обратите внимание, что для удобства я позволю x_array быть индексами первого измерения y_array, но в реальном анализе x_array наверняка можно прочитать извне.

Код

def linregress_3D(y_array):
    # y_array is a 3-D array formatted like (time,lon,lat)
    # The purpose of this function is to do linear regression using time series of data over each (lon,lat) grid box with consideration of ignoring np.nan
    # Construct x_array indicating time indexes of y_array, namely the independent variable.
    x_array=np.empty(y_array.shape)
    for i in range(y_array.shape[0]): x_array[i,:,:]=i+1 # This would be fine if time series is not too long. Or we can use i+yr (e.g. 2019).
    x_array[np.isnan(y_array)]=np.nan
    # Compute the number of non-nan over each (lon,lat) grid box.
    n=np.sum(~np.isnan(x_array),axis=0)
    # Compute mean and standard deviation of time series of x_array and y_array over each (lon,lat) grid box.
    x_mean=np.nanmean(x_array,axis=0)
    y_mean=np.nanmean(y_array,axis=0)
    x_std=np.nanstd(x_array,axis=0)
    y_std=np.nanstd(y_array,axis=0)
    # Compute co-variance between time series of x_array and y_array over each (lon,lat) grid box.
    cov=np.nansum((x_array-x_mean)*(y_array-y_mean),axis=0)/n
    # Compute correlation coefficients between time series of x_array and y_array over each (lon,lat) grid box.
    cor=cov/(x_std*y_std)
    # Compute slope between time series of x_array and y_array over each (lon,lat) grid box.
    slope=cov/(x_std**2)
    # Compute intercept between time series of x_array and y_array over each (lon,lat) grid box.
    intercept=y_mean-x_mean*slope
    # Compute tstats, stderr, and p_val between time series of x_array and y_array over each (lon,lat) grid box.
    tstats=cor*np.sqrt(n-2)/np.sqrt(1-cor**2)
    stderr=slope/tstats
    from scipy.stats import t
    p_val=t.sf(tstats,n-2)*2
    # Compute r_square and rmse between time series of x_array and y_array over each (lon,lat) grid box.
    # r_square also equals to cor**2 in 1-variable lineare regression analysis, which can be used for checking.
    r_square=np.nansum((slope*x_array+intercept-y_mean)**2,axis=0)/np.nansum((y_array-y_mean)**2,axis=0)
    rmse=np.sqrt(np.nansum((y_array-slope*x_array-intercept)**2,axis=0)/n)
    # Do further filteration if needed (e.g. We stipulate at least 3 data records are needed to do regression analysis) and return values
    n=n*1.0 # convert n from integer to float to enable later use of np.nan
    n[n<3]=np.nan
    slope[np.isnan(n)]=np.nan
    intercept[np.isnan(n)]=np.nan
    p_val[np.isnan(n)]=np.nan
    r_square[np.isnan(n)]=np.nan
    rmse[np.isnan(n)]=np.nan
    return n,slope,intercept,p_val,r_square,rmse

Пример вывода

Я использовал эту программу для тестирования двух трехмерных массивов с разрешением 227x3601x6301 пикселей, и она завершила работу в течение 20 минут, каждый менее чем за 10 минут.

Answer 2

На уровне numpy вы можете использовать np.vectorize .

Сначала определите сложную часть для каждого набора данных:

def compute(x,y):
        mask=~np.isnan(y)
        return linregress(x[mask],y[mask])

Затем определите функцию, которая будет обрабатывать ваши данные:

comp = np.vectorize(compute,signature="(k),(k)->(),(),(),(),()")

Затем примените, реорганизуя данные в соответствии с правилами вещания:

res = comp(x_array,rollaxis(y_array,0,3))

Наконец,

final=np.dstack(res) 

Теперь final[i,j] содержит пять параметров, возвращаемых linregress для пикселя (i,j).

Это примерно эквивалентно ответу метода панд, но в 2,5 раза быстрее.
Для решения проблемы 300x (изображение 100x100) требуется около 5 секунд, так что считайте один час для себя. Я не думаю, что это легко сделать лучше, так как время тратится в основном на функцию linregress.

Answer 3

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

https://hrishichandanpurkar.blogspot.com/2017/09/vectorized-functions-for-correlation.html

Я внес одно незначительное изменение (первая строка после #3), чтобы функция правильно учитывала различное количество значений 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

Answer 4

Я не уверен, как это увеличится (возможно, вы могли бы использовать dask ), но вот довольно простой способ сделать это с пандами DataFrame с помощью apply метод:

import pandas as pd
import numpy as np
from scipy.stats import linregress

# Independent variable: four time-steps of 1-dimensional data 
x_array = np.array([0.5, 0.2, 0.4, 0.4])

# Dependent variable: four time-steps of 3x3 spatial data
y_array = np.array([[[-0.2,   -0.2,   -0.3],
                     [-0.3,   -0.2,   -0.3],
                     [-0.3,   -0.4,   -0.4]],

                    [[-0.2,   -0.2,   -0.4],
                     [-0.3,   np.nan, -0.3],
                     [-0.3,   -0.3,   -0.4]],

                    [[np.nan, np.nan, -0.3],
                     [-0.2,   -0.3,   -0.7],
                     [-0.3,   -0.3,   -0.3]],

                    [[-0.1,   -0.3,   np.nan],
                     [-0.2,   -0.3,   np.nan],
                     [-0.1,   np.nan, np.nan]]])

def lin_regress(col):
    "Mask nulls and apply stats.linregress"
    col = col.loc[~pd.isnull(col)]
    return linregress(col.index.tolist(), col)

# Build the DataFrame (each index represents a pixel)
df = pd.DataFrame(y_array.reshape(len(y_array), -1), index=x_array.tolist())

# Apply a our custom linregress wrapper to each function, split the tuple into separate columns
final_df = df.apply(lin_regress).apply(pd.Series)

# Name the index and columns to make this easier to read
final_df.columns, final_df.index.name = 'slope, intercept, r_value, p_value, std_err'.split(', '), 'pixel_number'

print(final_df)

Выход:

                 slope  intercept   r_value       p_value   std_err
pixel_number                                                       
0             0.071429  -0.192857  0.188982  8.789623e-01  0.371154
1            -0.071429  -0.207143 -0.188982  8.789623e-01  0.371154
2             0.357143  -0.464286  0.944911  2.122956e-01  0.123718
3             0.105263  -0.289474  0.229416  7.705843e-01  0.315789
4             1.000000  -0.700000  1.000000  9.003163e-11  0.000000
5            -0.285714  -0.328571 -0.188982  8.789623e-01  1.484615
6             0.105263  -0.289474  0.132453  8.675468e-01  0.557000
7            -0.285714  -0.228571 -0.755929  4.543711e-01  0.247436
8             0.071429  -0.392857  0.188982  8.789623e-01  0.371154