我有一个矩阵M(比如r行和c列),我想根据它的邻居得到每个矩阵元素的“加权”和,并创建一个新的矩阵M2。邻居这个词可以在1的半径范围内(通常称为元胞自动机理论中的摩尔邻域),或者半径可以是1,例如2,3等。
对于矩阵M中的特定单元格,请说出中间的某个位置。让我们说位置(i,j);然后第(i,j)个小区有“八个”邻居给出,
(i-1, j-1), (i-1, j), (i-1, j+1), (i, j-1), (i, j+1), (i+1, j-1), (i+1, j), (i+1, j+1).
我想创建一个矩阵M2,它计算第(i,j)个单元加上其八个邻居的“加权”和。基于细胞之间的欧几里德距离进行加权。所以,例如,
exp(-sqrt(2))*M[i-1,j-1] + exp(-1)*M[i-1,j] + exp(-sqrt(2))*M[i-1,j+1] + exp(-1)*M[i,j-1] + M[i,j] + exp(-1)*M[i,j+1] + exp(-sqrt(2))*M[i+1,j-1] + exp(-1)*M[i+1,j] + exp(-sqrt(2))*M[i+1,j+1]
对于所有细胞重复相同的想法(沿着边界的细胞需要被特别处理,因为它们不一定具有八个相邻细胞)。上面的想法适用于半径1,但我想要开发的代码需要对任何半径都是通用的。
r <- 4
c <- 4
n <- r*c
(M <- matrix(1:n, r, c))
addresses <- expand.grid(x = 1:r, y = 1:c)
#Got this code in the same forum
z <- rbind(c(-1,0,1,-1,1,-1,0,1),c(-1,-1,-1,0,0,1,1,1))
get.neighbors <- function(rw) {
# Convert to absolute addresses
z2 <- t(z + unlist(rw))
# Choose those with indices within mat
b.good <- rowSums(z2 > 0)==2 & z2[,1] <= nrow(M) & z2[,2] <= ncol(M)
M[z2[b.good,]]
}
apply(addresses,1 , get.neighbors) # Returns a list with neighbors
M
本质上,半径= 1的M2必须是当前单元加上邻居的“加权”和。当前当前单元格的权重始终为1。
M = [ 1 5 9 13
2 6 10 14
3 7 11 15
4 8 12 16]
M2 = [ 5.033 13.803 .... ....
.... .... .... ....
.... .... .... ....
.... .... .... ....]
如何在R中获取矩阵M2?如果半径超过1怎么样?我想在两个for循环中发生加权,所以我可以在代码中进一步使用[i,j]单元格的计算加权和来关闭两个for循环。
我认为以下是你想要的加权和。我会以类似的方式找到邻居@ r2evans。
wtd_nbrs_sum <- function(input_matrix,
radius,
weight_matrix)
{
temp_1 <- matrix(data = 0,
nrow = nrow(x = input_matrix),
ncol = radius)
temp_2 <- matrix(data = 0,
nrow = radius,
ncol = ((2 * radius) + ncol(x = input_matrix)))
input_matrix_modified <- rbind(temp_2,
cbind(temp_1, input_matrix, temp_1),
temp_2)
output_matrix <- matrix(nrow = nrow(x = input_matrix),
ncol = ncol(x = input_matrix))
for(i in seq_len(length.out = nrow(x = input_matrix)))
{
for(j in seq_len(length.out = nrow(x = input_matrix)))
{
row_min <- (radius + (i - radius))
row_max <- (radius + (i + radius))
column_min <- (radius + (j - radius))
column_max <- (radius + (j + radius))
neighbours <- input_matrix_modified[(row_min:row_max), (column_min:column_max)]
weighted_sum <- sum(neighbours * weight_matrix)
output_matrix[i, j] <- weighted_sum
}
}
return(output_matrix)
}
r <- 4
c <- 4
n <- r*c
M <- matrix(data = 1:n,
nrow = r,
ncol = c)
R <- 1
wts <- matrix(data = c(exp(x = -sqrt(x = 2)), exp(x = -1), exp(x = -sqrt(x = 2)), exp(x = -1), 1, exp(x = -1), exp(x = -sqrt(x = 2)), exp(x = -1), exp(x = -sqrt(x = 2))),
nrow = 3,
ncol = 3)
wtd_nbrs_sum(input_matrix = M,
radius = R,
weight_matrix = wts)
#> [,1] [,2] [,3] [,4]
#> [1,] 5.033856 13.80347 24.16296 23.89239
#> [2,] 8.596195 20.66391 34.43985 32.84175
#> [3,] 11.186067 24.10789 37.88383 35.43163
#> [4,] 9.748491 19.86486 30.22435 28.60703
由reprex package创建于2019-03-24(v0.2.1)
希望这可以帮助。
编辑包括加权和。
可能有一个非常巧妙的技巧,但最直接(和可维护)的方式可能是一个简单的两个for
循环实现。
M1 <- matrix(1:16, nr=4)
M1
# [,1] [,2] [,3] [,4]
# [1,] 1 5 9 13
# [2,] 2 6 10 14
# [3,] 3 7 11 15
# [4,] 4 8 12 16
代码:
get_neighbors <- function(M, radius = 1) {
M2 <- M
M2[] <- 0
nr <- nrow(M)
nc <- ncol(M)
eg <- expand.grid((-radius):radius, (-radius):radius)
eg$wt <- exp(-sqrt(abs(eg[,1]) + abs(eg[,2])))
for (R in seq_len(nr)) {
for (C in seq_len(nc)) {
ind <- cbind(R + eg[,1], C + eg[,2], eg[,3])
ind <- ind[ 0 < ind[,1] & ind[,1] <= nr &
0 < ind[,2] & ind[,2] <= nc,, drop = FALSE ]
M2[R,C] <- sum(M[ind[,1:2, drop=FALSE]] * ind[,3])
}
}
M2
}
get_neighbors(M1, 1)
# [,1] [,2] [,3] [,4]
# [1,] 5.033856 13.80347 24.16296 23.89239
# [2,] 8.596195 20.66391 34.43985 32.84175
# [3,] 11.186067 24.10789 37.88383 35.43163
# [4,] 9.748491 19.86486 30.22435 28.60703
同样的事情,半径为2:
get_neighbors(M1, 2)
# [,1] [,2] [,3] [,4]
# [1,] 12.44761 25.64963 31.73247 32.70974
# [2,] 18.57765 35.96237 43.33862 43.51911
# [3,] 20.09836 37.80643 45.18268 45.03982
# [4,] 17.51314 31.88500 37.96784 37.77527
一个简单的测试,如果使用半径0,那么M1和M2应该是相同的(它们是)。
注意:这通常在基数R中表现不错,没有花哨使用apply
或其表兄弟。由于这是一个非常简单的启发式算法,因此可以很容易地使用Rcpp实现更快。