4

私は、3 つの空間次元と時間次元 (x、y、z、t) で構成される 4 次元配列で構成される R のデータに取り組んでいます。いくつかの分析では、一連の空間座標 x、y、z の時間次元のすべてのデータを取得したいと考えています。ここまでは、 which 関数を使用して、対象の空間位置のインデックスを取得しました。しかし、空間位置に対応する時間次元ですべての関連データを取得しようとすると、洗練された R ソリューションを見つけることができず、移植された MATLAB 関数である repmat を使用することになりました。

a4d <- array(rnorm(10000), rep(10,4)) #x, y, z, t

#arbitrary set of 3d spatial indices x, y, z (here, using high values at first timepoint)
indices <- which(a4d[,,,1] > 2, arr.ind=TRUE)
str(indices)

# int [1:20, 1:3] 10 2 6 5 8 2 6 8 2 10 ...
# - attr(*, "dimnames")=List of 2
# ..$ : NULL
# ..$ : chr [1:3] "dim1" "dim2" "dim3"

#Now, I would like to use these indices to get data x, y, z for all t

#Intuitive, but invalid, syntax (also not clear what the structure of the data would be)
#a4d[indices,]

#Ugly, but working, syntax
library(pracma)

#number of timepoints
nt <- dim(a4d)[4]

#create a 4d lookup matrix
lookup <- cbind(repmat(indices, nt, 1), rep(1:nt, each=nrow(indices)))

#obtain values at each timepoint for indices x, y, z
result <- cbind(lookup, a4d[lookup])

このソリューションは、指定された目的には問題なく機能しますが、概念的には醜いようです。理想的には、最後に 2 次元の行列 (インデックス x 時間) が必要です。したがって、この場合、ルックアップに 20 個の x、y、z 座標、および 10 個のタイムポイントがある場合、行がインデックスの各行を表す 20 x 10 行列が理想的です (x、y、z を保持する必要はありません)。 、値は必ず)、各列は時点です。

Rでこれを行う良い方法はありますか? do.call("[", list ... などをいじり、outer と prod を使用しましたが、期待どおりに機能しませんでした。

ご提案ありがとうございます。マイケル

4

3 に答える 3

7

私はあなたが探していると思います：

apply(a4d, 4, `[`, indices)

そして、結果が一致することを確認するには:

result1 <- matrix(result[,5], ncol = 10)
result2 <- apply(a4d, 4, `[`, indices)
identical(result1, result2)
# [1] TRUE
于 2012-07-17T02:21:36.190 に答える
1

私はおそらく何かが欠けていますが、あなたはただ欲しくないですa4d[indices[,1],indices[,2],indices[,3],]か？

于 2012-07-17T01:55:58.337 に答える
1

各次元で個別にアクセスすると、@tilo-wiklund または期待どおりに機能しません。10 タイム ステップで 23 行ではなく、結果は 10 タイム ステップで 23x23x23 の立方体になります。

r.idvdim <- a4d[indices[,1],indices[,2],indices[,3],]
r.apply  <- apply(a4d, 4, `[`, indices)
r.cbind  <- matrix(a4d[lookup],ncol=nt)

dim(r.idvdim)     # [1] 23 23 23 10
dim(r.apply)      # [1] 23 10
dim(r.cbind)      # [1] 23 10
于 2012-09-17T13:26:55.383 に答える