Python中如何获取矩阵中带有缺失值的邻居节点的平均值?

3

我有一个非常大的矩阵,因此不想通过遍历每一行和每一列来求和。

a = [[1,2,3],[3,4,5],[5,6,7]]
def neighbors(i,j,a):
    return [a[i][j-1], a[i][(j+1)%len(a[0])], a[i-1][j], a[(i+1)%len(a)][j]]
[[np.mean(neighbors(i,j,a)) for j in range(len(a[0]))] for i in range(len(a))]

这段代码对于3x3或较小范围的矩阵运行良好,但对于像2k x 2k这样的大型矩阵来说不可行。并且,如果矩阵中的任何值缺失或为na等情况时,此代码也无法工作。如果任何相邻值为na,则跳过该相邻值以获取平均值。

3个回答

6

第一步

假设您希望在一个输入数组中获得窗口大小为3 x 3,并仅考虑北、西、东、南四个相邻元素的滑动窗口平均值。

对于这种情况,可以使用signal.convolve2d函数和适当的卷积核。最后,您需要将这些总和除以卷积核中的1的数量,即kernel.sum(),因为只有这些元素对总和产生了贡献。以下是实现方法-

import numpy as np
from scipy import signal

# Inputs
a = [[1,2,3],[3,4,5],[5,6,7],[4,8,9]]

# Convert to numpy array
arr = np.asarray(a,float)    

# Define kernel for convolution                                         
kernel = np.array([[0,1,0],
                   [1,0,1],
                   [0,1,0]]) 

# Perform 2D convolution with input data and kernel 
out = signal.convolve2d(arr, kernel, boundary='wrap', mode='same')/kernel.sum()

第二步

这个步骤和第一步的假设相同,只是我们希望在仅有零元素的邻域中寻找平均值,并用这些平均值替换它们。

方法 #1:以下是使用手动选择卷积方法完成此操作的一种方式 -

import numpy as np

# Convert to numpy array
arr = np.asarray(a,float)    

# Pad around the input array to take care of boundary conditions
arr_pad = np.lib.pad(arr, (1,1), 'wrap')

R,C = np.where(arr==0)   # Row, column indices for zero elements in input array
N = arr_pad.shape[1]     # Number of rows in input array

offset = np.array([-N, -1, 1, N])
idx = np.ravel_multi_index((R+1,C+1),arr_pad.shape)[:,None] + offset

arr_out = arr.copy()
arr_out[R,C] = arr_pad.ravel()[idx].sum(1)/4

样例输入,输出 -

In [587]: arr
Out[587]: 
array([[ 4.,  0.,  3.,  3.,  3.,  1.,  3.],
       [ 2.,  4.,  0.,  0.,  4.,  2.,  1.],
       [ 0.,  1.,  1.,  0.,  1.,  4.,  3.],
       [ 0.,  3.,  0.,  2.,  3.,  0.,  1.]])

In [588]: arr_out
Out[588]: 
array([[ 4.  ,  3.5 ,  3.  ,  3.  ,  3.  ,  1.  ,  3.  ],
       [ 2.  ,  4.  ,  2.  ,  1.75,  4.  ,  2.  ,  1.  ],
       [ 1.5 ,  1.  ,  1.  ,  1.  ,  1.  ,  4.  ,  3.  ],
       [ 2.  ,  3.  ,  2.25,  2.  ,  3.  ,  2.25,  1.  ]])

为了处理边界条件,填充有其他选项。请查看numpy.pad以获取更多信息。

方法 #2:这是之前在Shot #1中列出的基于卷积的方法的修改版。这与之前的方法相同,只是在最后,我们选择性地用卷积输出替换零元素。下面是代码 -

import numpy as np
from scipy import signal

# Inputs
a = [[1,2,3],[3,4,5],[5,6,7],[4,8,9]]

# Convert to numpy array
arr = np.asarray(a,float)

# Define kernel for convolution                                         
kernel = np.array([[0,1,0],
                   [1,0,1],
                   [0,1,0]]) 

# Perform 2D convolution with input data and kernel 
conv_out = signal.convolve2d(arr, kernel, boundary='wrap', mode='same')/kernel.sum()

# Initialize output array as a copy of input array
arr_out = arr.copy()

# Setup a mask of zero elements in input array and 
# replace those in output array with the convolution output
mask = arr==0
arr_out[mask] = conv_out[mask]

备注:如果输入数组中零元素较少,则首选方法1,否则请使用方法2


能否仅针对矩阵中特定点获取平均值,例如值为“8”的元素? - Aamirkhan
1
@Aamirkhan 你的意思是对于一个特定的坐标列表/数组,将ij作为输入,然后我们会计算出这些坐标周围邻近区域的平均值? - Divakar
抱歉给您带来不便。实际上,矩阵中的0值元素被识别为缺失值。因此,我想通过获取其邻域元素的平均值来填充它。所以我想要获取仅值为0的元素的邻居的平均值。 - Aamirkhan
@paddyg 当然可以,那是一种方法。但是,如果输入数组中的零元素非常少,那么这种方法效率不高。因此,我在解决方案中实现了一种“手动卷积”方法,即Shot#2。 - Divakar
@Divakar,虽然我投票支持了你的答案,但我在下面发布了一些有趣的timeit()结果... - paddyg
显示剩余5条评论

3

这是对@Divakar答案下评论的补充(而不是一个独立的答案)。

出于好奇,我尝试了不同的“伪”卷积方法来与scipy卷积进行比较。最快的是%(模数)包裹方法,这让我感到惊讶:显然numpy在其索引方面做了一些聪明的事情,但显然不需要填充将节省时间。

fn3 -> 9.5毫秒,fn1 -> 21毫秒,fn2 -> 232毫秒

import timeit

setup = """
import numpy as np
from scipy import signal
N = 1000
M = 750
P = 5 # i.e. small number -> bigger proportion of zeros
a = np.random.randint(0, P, M * N).reshape(M, N)
arr = np.asarray(a,float)"""

fn1 = """ 
arr_pad = np.lib.pad(arr, (1,1), 'wrap')
R,C = np.where(arr==0)
N = arr_pad.shape[1]
offset = np.array([-N, -1, 1, N])
idx = np.ravel_multi_index((R+1,C+1),arr_pad.shape)[:,None] + offset
arr[R,C] = arr_pad.ravel()[idx].sum(1)/4"""

fn2 = """
kernel = np.array([[0,1,0],
                   [1,0,1],
                   [0,1,0]]) 
conv_out = signal.convolve2d(arr, kernel, boundary='wrap', mode='same')/kernel.sum()
mask = arr == 0.0
arr[mask] = conv_out[mask]"""

fn3 = """ 
R,C = np.where(arr == 0.0)
arr[R, C] = (arr[(R-1)%M,C] + arr[R,(C-1)%N] + arr[R,(C+1)%N] + arr[(R+1)%M,C]) / 4.0
"""

print(timeit.timeit(fn1, setup, number = 100))
print(timeit.timeit(fn2, setup, number = 100))
print(timeit.timeit(fn3, setup, number = 100))

需要填充来处理边界条件,然后使用填充数组中的元素。因此,这些肯定会增加开销和通用情况下的运行时间。但是,原地替换是很好的工作! - Divakar

2
使用 numpyscipy.ndimage,您可以应用一个“足印”,该足印定义了您查找每个元素的邻居并将函数应用于这些邻居的位置。
import numpy as np
import scipy.ndimage as ndimage

# Getting neighbours horizontally and vertically,
#   not diagonally
footprint = np.array([[0,1,0],
                      [1,0,1],
                      [0,1,0]])
a = [[1,2,3],[3,4,5],[5,6,7]]
# Need to make sure that dtype is float or the
#   mean won't be calculated correctly
a_array = np.array(a, dtype=float)

# Can specify that you want neighbour selection to
#   wrap around at the borders
ndimage.generic_filter(a_array, np.mean, 
                       footprint=footprint, mode='wrap')
Out[36]: 
array([[ 3.25,  3.5 ,  3.75],
       [ 3.75,  4.  ,  4.25],
       [ 4.25,  4.5 ,  4.75]])

能否仅针对矩阵中特定点获取平均值,例如值为“8”的元素? - Aamirkhan

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接