Gauss矩阵迭代

Question

假设我有一个矩阵：

> mat <- matrix(c( rep(1:12,1)), 4, 4)
> mat
     [,1] [,2] [,3] [,4]
[1,]    1    5    9    1
[2,]    2    6   10    2
[3,]    3    7   11    3
[4,]    4    8   12    4

我希望使用被包围元素的平均值来替换中间的4个值；我可以使用以下方法来实现：

for (i in 1:1000){
  mat[2,2]=  (mat[1,2] + mat[2,1] + mat[3,2] + mat[2,3])/4
  mat[2,3]=  (mat[2,2] + mat[2,4] + mat[1,3] + mat[3,3])/4
  mat[3,2]=  (mat[3,1] + mat[3,3] + mat[2,2] + mat[4,2])/4
  mat[3,3]=  (mat[3,2] + mat[3,4] + mat[4,3] + mat[2,3])/4
  print(mat)
}

我的问题是:我如何将它转化为一个函数，以便我可以直接将它应用于选定的矩阵元素(即mat[2:3,2:3])？

Answer 1

GSfun <- function(M, r, c) {sum(matrix(c(0,1,0,
                                         1,0,1,
                                         0,1,0), 3) * 
                                 M[r+(-1:1),c+(-1:1)]) /4 }

for (i in 2:3){ for (j in 2:3){  mat[i,j] <- GSfun(mat, i,j)}}
mat
     [,1] [,2]  [,3] [,4]
[1,]    1    5  9.00    1
[2,]    2    6  7.00    2
[3,]    3    7  7.25    3
[4,]    4    8 12.00    4
for (i in 2:3){ for (j in 2:3){  mat[i,j] <- GSfun(mat, i,j)}}
mat
     [,1]  [,2]    [,3] [,4]
[1,]    1 5.000  9.0000    1
[2,]    2 5.250  5.8750    2
[3,]    3 5.875  6.6875    3
[4,]    4 8.000 12.0000    4

迭代应用，直到绝对误差之和小于0.01为止：

 maxiter = 20; n=1; repeat { n=n+1; 
     diverg <- sum(abs(mat[2:3, 2:3])); 
     for (i in 2:3){ for (j in 2:3){  mat[i,j] <- GSfun(mat, i,j)}}; 
     if ( n>maxiter | (diverg - sum(abs(mat[2:3, 2:3])) < 0.01) ) {break}  }
 n
# [1] 8
 mat
#------------
     [,1]     [,2]      [,3] [,4]
[1,]    1 5.000000  9.000000    1
[2,]    2 4.500732  5.500366    2
[3,]    3 5.500366  6.500183    3
[4,]    4 8.000000 12.000000    4

Answer 2

试试这个：

# a matrix with slightly 'simpler' numbers
set.seed(1)
mm <- matrix(sample(1:3, 16, replace = TRUE), ncol = 4)
mm

#      [,1] [,2] [,3] [,4]
# [1,]    1    1    2    3
# [2,]    2    3    1    2
# [3,]    2    3    1    3
# [4,]    3    2    1    2

f_replace <- function(mm, rows, cols){
  mm2 <- mm
  for(i in rows){
    for(j in cols){
      rows_around <- i + c(0, 0, -1, 1)
      cols_around <- j + c(-1, 1, 0, 0)
      id <- cbind(rows_around, cols_around)
      mm2[i, j] <- mean(mm[id])
    }
  }
  return(mm2)
}

f_replace(mm, 2:3, 2:3)
#      [,1] [,2] [,3] [,4]
# [1,]    1 1.00    2    3
# [2,]    2 1.75    2    2
# [3,]    2 2.00    2    3
# [4,]    3 2.00    1    2

# a loop where you can plug in your n of choice
for(i in n){
  mm <- f_replace(mm, 2:3, 2:3)
}