RでMATLABのrepmat実装
—– 2010/12/25 追記 ————–
コードを修正したんで新しくエントリー書きました
———————————–
個人的にMATLABよりも断然R派の私ですが，MATLABをちょっとだけ使ってみて，デフォルトで用意されている関数で唯一便利だと思う関数があります．
repmat
repmat(A,m,n) は行列 A を行方向に m 回，列方向に n 回繰り返した行列を返します．
例えばこんな感じです．
>> A = [1:3; 4:6]
A =
1 2 3
4 5 6
>> repmat(A, 2, 3)
ans =
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
repmat(A, [m1, m2, …, mN]) にすると，repmat(A, m1, m2) で作成した行列を各次元に指定した回数だけ繰り返した多次元配列を返します．
例えばこんな感じです．
>> repmat(A, [2 3 2])
ans(:,:,1) =
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
ans(:,:,2) =
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
Aが3次元以上の配列の場合も対応しているようです．
>> repmat(repmat(A, [2 3 2]), 2, 1)
ans(:,:,1) =
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
ans(:,:,2) =
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
Rでもこういった関数があればたまには重宝しそうです．
そんなわけで実装してみました．kroneckerを使う方法よりは速いのと，中途半端にAが3次元以上の配列の場合にも対応しているのが特徴です．
以下の実装ではAが3次元以上であり，第2引数も3次元以上の場合は非対応です．けっこう頑張りましたがどうすれば互換性を持たせて対応させられるかわかりません・・・．しかもMATLABの公式サイトのrepmatの説明を見る限り，実装上3次元以上にも対応しているだけであり，3次元以上の配列は想定していない気がします．
repmat <- function(A, m, n = m) {
if(missing(m))
stop("Requires at least 2 inputs.")
if(is.vector(A))
A <- t(matrix(A))
d <- dim(A)
nr <- d[1]
nc <- d[2]
if(length(m) <= 2) {
m[2] <- ifelse(is.na(m[2]), n, m[2])
if(length(d) > 2) {
## 2次元になるまで再帰で次元を削減
ret <- array(dim = c(nr * m[1], nc * m[2], d[-(1:2)]))
for(i in 1:d[length(d)]){
sep <- paste(rep(",", length(d) - 1), collapse="")
eval(parse(text = sprintf("ret[%si] <- repmat(A[%si], m[1], m[2])", sep, sep)))
}
} else {
tmpA <- matrix(rep(t(A), m[1]), nrow = nr * m[1], byrow = TRUE)
return(matrix(rep(tmpA, m[2]), nrow = nr * m[1], ncol = nc * m[2]))
}
ret
} else {
if(length(d) > 2)
stop("Doesn't support these arguments.") # MATLABではサポートされている
tmpA <- matrix(rep(t(A), m[1]), nrow = nr * m[1], byrow = TRUE)
tmpA <- matrix(rep(tmpA, m[2]), nrow = nr * m[1], ncol = nc * m[2])
array(rep(tmpA, prod(m[-(1:2)])), c(nr * m[1], nc * m[2], m[-(1:2)]))
}
}
実際に使ってみます．
> (A <- matrix(1:6, 2, byrow=T))
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
> repmat(A, 2, 3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 1 2 3 1 2 3 1 2 3
[2,] 4 5 6 4 5 6 4 5 6
[3,] 1 2 3 1 2 3 1 2 3
[4,] 4 5 6 4 5 6 4 5 6
> repmat(A, c(2, 3, 2))
, , 1
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 1 2 3 1 2 3 1 2 3
[2,] 4 5 6 4 5 6 4 5 6
[3,] 1 2 3 1 2 3 1 2 3
[4,] 4 5 6 4 5 6 4 5 6
, , 2
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 1 2 3 1 2 3 1 2 3
[2,] 4 5 6 4 5 6 4 5 6
[3,] 1 2 3 1 2 3 1 2 3
[4,] 4 5 6 4 5 6 4 5 6
> repmat(repmat(A, c(2, 3, 2)), 2, 1)
, , 1
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 1 2 3 1 2 3 1 2 3
[2,] 4 5 6 4 5 6 4 5 6
[3,] 1 2 3 1 2 3 1 2 3
[4,] 4 5 6 4 5 6 4 5 6
[5,] 1 2 3 1 2 3 1 2 3
[6,] 4 5 6 4 5 6 4 5 6
[7,] 1 2 3 1 2 3 1 2 3
[8,] 4 5 6 4 5 6 4 5 6
, , 2
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 1 2 3 1 2 3 1 2 3
[2,] 4 5 6 4 5 6 4 5 6
[3,] 1 2 3 1 2 3 1 2 3
[4,] 4 5 6 4 5 6 4 5 6
[5,] 1 2 3 1 2 3 1 2 3
[6,] 4 5 6 4 5 6 4 5 6
[7,] 1 2 3 1 2 3 1 2 3
[8,] 4 5 6 4 5 6 4 5 6
第2引数が3次元以上の場合にも対応させようとして非常に不毛な時間を過ごした気がする…