我有一个独立变量值的一维数组 (x_array)，它与具有多个时间步长的 3D numpy 空间数据数组 (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]]])

我想计算每个像素的线性回归并获得 y_array 中每个 xy 像素的 R 平方、P 值、截距和斜率，其值为x_array 中的每个时间步作为我的自变量。

我可以 reshape 以获取表单中的数据，以将其输入到矢量化且快速的 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 值、截距和斜率)。

answer here使用 scipy.stats import linregress 计算我需要的统计数据，并建议通过屏蔽这些 NaN 值来避免 NaN 问题。然而，这个例子是针对两个一维数组的，我不知道如何对我的情况应用类似的掩码方法，在这种情况下，y_array_reshaped 中的每一列都有一组不同的 NaN 值。

我的问题:如何计算一个包含许多 NaN 值的大型多维数组 (300 x 3000 x 3000) 中每个像素的回归统计数据，速度相当快，矢量化方式？

期望的结果:y_array 中每个像素的回归统计值(例如 R 平方)的 3 x 3 数组，即使该像素包含 NaN 时间序列中某个点的值

最佳答案

上面评论中提到的这篇博文包含一个非常快速的向量化函数，用于 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

关于python - 如何对包含 NaN 的大型多维数组中的每个像素应用线性回归？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/52108417/