R_winBUGS_structure のバックアップ(No.9)

バックアップ一覧
差分を表示
現在との差分を表示
ソースを表示
R_winBUGS_structure へ行く。
- 1 (2008-09-24 (水) 16:30:48)
- 2 (2008-09-24 (水) 17:49:27)
- 3 (2008-09-25 (木) 08:18:16)
- 4 (2008-09-25 (木) 13:30:10)
- 5 (2008-09-26 (金) 12:18:13)
- 6 (2008-09-26 (金) 22:19:34)
- 7 (2008-09-29 (月) 11:37:22)
- 8 (2008-09-30 (火) 16:02:10)
- 9 (2008-10-07 (火) 21:25:11)
- 10 (2008-10-09 (木) 18:24:20)

以下，WinBUGSデータの配列を R 上の JAGS で扱えるようにするための検討をしていたのだが，coda パッケージにすでに実装されていることに気がついた．すなわちbugs2jags()関数である．

bugs2jags("winbugsData.txt", "jagsData.txt")

と使えばよろしい．．．

以下，消すのも何なので残すが，参照の必要なし．

3次元配列
4次元配列

structure 関数が R では，縦方向に行列を埋めていくのに対して，winBUGS では横方向にデータを埋めていくのであった．例えば以下は

structure(
   .Data = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
   .Dim=c(5,2)

winBUGS では

[1] 1 2 3 4 5
[2] 6 7 8 9 10

となるが(winBUGSのメニュー[Info]-[Node Info]-[y]-[values]で確認してみた)，R では

[1] 1 3 5 7 9 
[2] 2 4 6 8 10

となる．

_ 3次元配列

X1 <- structure(.Data = 1:24, .Dim=c(2,3,4))

はWinBUGSでは,最後に指定された3次元の添字1:4を繰り返し，次に最後から二つ目，つまり2次元の添字1:3を繰り返して埋めていく．

y[1,1,1]      1.0
y[1,1,2]      2.0
y[1,1,3]      3.0
y[1,1,4]      4.0
y[1,2,1]      5.0
y[1,2,2]      6.0
y[1,2,3]      7.0
...
y[2,3,1]      21.0
y[2,3,2]      22.0
y[2,3,3]      23.0
y[2,3,4]      24.0

となる．これはRの配列で言うと

> X2
, , 1
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]   13   14   15
, , 2
     [,1] [,2] [,3]
[1,]    4    5    6
[2,]   16   17   18
, , 3
     [,1] [,2] [,3]
[1,]    7    8    9
[2,]   19   20   21
, , 4
     [,1] [,2] [,3]
[1,]   10   11   12
[2,]   22   23   24

ということになるが，先のstrucute関数をそのままRで実行すると

, , 1
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
, , 2 
     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12
, , 3
     [,1] [,2] [,3]
[1,]   13   15   17
[2,]   14   16   18
, , 4
     [,1] [,2] [,3]
[1,]   19   21   23
[2,]   20   22   24

となる．winBUGSのstructure()関数と同じ配列を作成するためには以下のようにしなければならない．

> X2 <- structure(.Data =1,13,2,14,3,15,4,16,5,17,6,18,7,19,8,
                               20,9,21,10,22,11,23,12,24, .Dim=c(2,3,4))

すなわち以下のベクトルを

> X1 <- 1:24
> X1
 1   2    3   4  5   6  7   8  9 10 11 12
 13 14 15 16 17 18 19 20 21 22 23 24
> Y1
 1 13  2 14  3 15  4 16  5 17  6 18  7 19  8 20  9 21 10 22
 11 23 12 24

のように並び替える必要がある．つまり3*4個ずらして並べ直すのである．

for(i in 1: length(X1)){
 Y1 <- c(Y1,  X1[seq(i, length(X1), 12) ])
 if(length(Y1) >= length(X1)) break
}

具体例は以下

y=structure(
 .Data = c(75, 57, 10, 66, 74, 40, 52, 60, 54,  9, 64, 58,
           63, 31, 63, 58, 67, 44, 20, 53, 65, 10, 42, 20), 
 .Dim = c(2, 3, 4)))

はwinBUGSでは,最後に指定された3次元の添字を繰り返し，次に最後から二つ目，つまり2次元の添字を加算して，埋めていく．

model compiled
y[1,1,1]      75.0
y[1,1,2]      57.0
y[1,1,3]      10.0
y[1,1,4]      66.0
y[1,2,1]      74.0
y[1,2,2]      40.0
y[1,2,3]      52.0
y[1,2,4]      60.0
y[1,3,1]      54.0
y[1,3,2]      9.0
y[1,3,3]      64.0
y[1,3,4]      58.0
y[2,1,1]      63.0
y[2,1,2]      31.0
y[2,1,3]      63.0
y[2,1,4]      58.0
y[2,2,1]      67.0
y[2,2,2]      44.0
y[2,2,3]      20.0
y[2,2,4]      53.0
y[2,3,1]      65.0
y[2,3,2]      10.0
y[2,3,3]      42.0
y[2,3,4]      20.0

と展開される。これはRの配列で書くと

ということになるが，しかしRで同じstructure()関数を実行すると

y=structure(
 .Data = c(75, 57, 10, 66, 74, 40, 52, 60, 54,  9, 64, 58,
           63, 31, 63, 58, 67, 44, 20, 53, 65, 10, 42, 20), 
 .Dim = c(2, 3, 4)))

結果は

, , 1
      [,1] [,2] [,3]
[1,]   75   10   74
[2,]   57   66   40
 , , 2
      [,1] [,2] [,3]
[1,]   52   54   64
[2,]   60    9   58
 , , 3
      [,1] [,2] [,3]
[1,]   63   63   67
[2,]   31   58   44
 , , 4
      [,1] [,2] [,3]
[1,]   20   65   42
[2,]   53   10   20

となってしまう．winBUGSのstructure()関数と同じ配列を作成するためには以下のようにしなければならない．

y = structure(
      .Data = c(75, 63, 57, 31, 10, 63, 66, 58, 74, 67, 40, 44, 52, 
                 20, 60, 53, 54, 66, 9, 10, 64, 42, 58, 20), 
       .Dim = c(2, 3, 4))

すなわち

75, 57, 10, 66, 74, 40, 52, 60, 54,  9, 64, 58, 
63, 31, 63, 58, 67, 44, 20, 53, 65, 10, 42, 20

を

75, 63, 57, 31, 10, 63, 66, 58, 74, 67, 40, 44, 
52, 20, 60, 53, 54, 66, 9, 10, 64, 42, 58, 20

と並び替える．3*4 個ずつずれるのである．だから，手っ取り早くは次のようにして並び替えができる．

y <- numeric()
for(i in 1: length(x)){
  y <- c(y,  x[seq(i, length(x), 12) ])
  if(length(y) >= length(x)) break
}

また先のwinBUGSの配列は R のデータフレーム流に書くと

> df.array <- data.frame(I = c(rep("1",12), rep("2",12)), 
             J = c( rep("1",4), rep("2", 4), rep("3", 4)), 
             K = c(75, 57, 10, 66, 74, 40, 52, 60, 54, 9, 
                     64, 58, 63, 31, 63, 58, 67, 44, 20, 53, 
                     65, 10, 42, 20))
> df.array
  I J  K
1  1 1 75
2  1 1 57
3  1 1 10
4  1 1 66
5  1 2 74
6  1 2 40
7  1 2 52
8  1 2 60
9  1 3 54
10 1 3  9
11 1 3 64
12 1 3 58
13 2 1 63
14 2 1 31
15 2 1 63
16 2 1 58
17 2 2 67
18 2 2 44
19 2 2 20
20 2 2 53
21 2 3 65
22 2 3 10
23 2 3 42
24 2 3 20

ということになる．

http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/manual14.pdf の17ページに詳しく書いてある．が，このファイルの説明は R ではうまくいかない．変換は，行列であれば転置すれば良いだけだが，配列だとやっかいだな．作成時にstructure()関数に引数を与えることで，横方向に埋めることは可能なのだろうか？

M <- matrix(1:10, ncol = 2)
str(M)
# int [1:5, 1:2] 1 2 3 4 5 6 7 8 9 10
str(t(M))
# int [1:2, 1:5] 1 6 2 7 3 8 4 9 5 10
dput(M, file = "M.txt")
# structure(1:10, .Dim = c(5L, 2L))
dput(t(M), file = "M.txt")
# structure(c(1L, 6L, 2L, 7L, 3L, 8L, 4L, 9L, 5L, 10L), 
  .Dim = c(2L,  5L))

ちなみに上の出力の数字の後にLがあるが，これは整数値として保存されていることを意味する．

Rのメーリングリストで [R] Bug in "is" ? という話題があった．その2008年9月25日の回答に以下のようなものがあった．

Mathematicians are concerned with properties of numbers, computer scientists are concerned with how numbers are stored (and statisticians when doing statistics are more concerned with data than numbers). R is an implementation of the S programming language (along with many tools written in that language) so fits in more with the computer scientist view than the mathematical view. So, is.integer is telling you about how 7 is stored, not the property of the number 7. If we write 7L then we tell R/S that we want 7 stored as an integer, if we write 7. or 7.0 then we tell R/S to store it as double precision, but with 7 it has to guess which we want, and since the real numbers are closed to more operations than the integers (and double precision is the chosen approximation to real), it seems the more practical default. Greg Snow

↑

_ 4次元配列

さて，

X=structure( 
.Data=c(1,2,3,4,5,6,7,8,9,10,
11,12,13,14,15,16,17,18,19,20,
21,22,23,24,25, 26,27,28,29,30,
31,32,33,34,35,36,37,38,39,40,
41,42,43,44,45,46,47,48,49,50,
51,52,53,54,55,56,57,58,59,60,
61,62,63,64,65,66,67,68,69,70,
71,72,73,74,75,76,77,78,79,80,
81,82,83,84,85,86,87,88,89,90,
91,92,93,94,95,96,97,98,99,100,
101,102,103,104,105,106,107,108,109,110,
111,112,113,114,115,116,117,118,119,120),
    .Dim=c(5,2,3,4)
 )

は WinBUGSではこんな風に展開される．指定された最後の4次元の添字を繰り返し，次に最後から二つ目，３次元の添字を繰り返す．

data loaded 
X[1,1,1,1]      1.0
X[1,1,1,2]      2.0
X[1,1,1,3]      3.0
X[1,1,1,4]      4.0
X[1,1,2,1]      5.0
X[1,1,2,2]      6.0
X[1,1,2,3]      7.0
X[1,1,2,4]      8.0
X[1,1,3,1]      9.0
X[1,1,3,2]      10.0
X[1,1,3,3]      11.0
X[1,1,3,4]      12.0
X[1,2,1,1]      13.0
X[1,2,1,2]      14.0
... 
X[5,2,1,1]      109.0
X[5,2,1,2]      110.0
X[5,2,1,3]      111.0
X[5,2,1,4]      112.0
X[5,2,2,1]      113.0
X[5,2,2,2]      114.0
X[5,2,2,3]      115.0
X[5,2,2,4]      116.0
X[5,2,3,1]      117.0
X[5,2,3,2]      118.0
X[5,2,3,3]      119.0
X[5,2,3,4]      120.0

で，これをRで実行すると,以下のように3時限目の添字を最初に繰り返して，次に4次元目の添字を繰り返している．

, , 1, 1
     [,1] [,2]
[1,]    1    6
[2,]    2    7
[3,]    3    8
[4,]    4    9
[5,]    5   10
, , 2, 1
     [,1] [,2]
[1,]   11   16
[2,]   12   17
[3,]   13   18
[4,]   14   19
[5,]   15   20
, , 3, 1
     [,1] [,2]
[1,]   21   26
[2,]   22   27
[3,]   23   28
[4,]   24   29
[5,]   25   30
...
, , 3, 4
     [,1] [,2]
[1,]  111  116
[2,]  112  117
[3,]  113  118
[4,]  114  119
[5,]  115  120

アールメカブ