Applied Statistics and Data Mining

マルコフ連鎖とマルコフ決定過程は、状態が多次元の場合にも応用できる。以下では、状態が２次元の場合を例に、１次元に変換して、マルコフ連鎖のパッケージmarkovchainとマルコフ決定過程のパッケージMDPtoolboxを用いる手順を説明する。〔現時点でこの原稿は、データの変換と状態遷移行列の推定については省略し、マルコフ決定過程の適用の準備を中心に述べている。〕

状態が2次元ベクトル(X.1, X.2)で表される場合を考えたい。 X.1は空間S.1に含まれるとする。

S.1 = c(0:5)

X.1は空間S.1に含まれるとする。

S.2 = c(0:5)

それぞれの状態空間を表示させてみる。

S.1
S.2

これらのすべての組み合わせの一覧を生成するには、関数expand.gridを用いる。

expand.grid(S.1,S.2)

これを直積状態空間という。 関数expand.gridは、直積集合(複数の集合のすべての要素の組み合わせ、からなる集合)を生成する。 3つ以上の集合でも動作する。

expand.grid(S.1, S.1, S.1, S.1)

が、ここでは2つで動いてくれれば十分である。

S.1とS.2は異なっていてもいいが、このページでは次の2つを仮定する。

• 0から始まること
• 一番大きな数字が状態数-1であること

この直積空間のすべての組み合わせには、一番左に表示される数字が、新しく付与される。 この数字を「状態」と思うと、直積空間を１次元空間とみなせる。

> expand.grid(c(0:4),c(0:4)) # 実行例
   Var1 Var2
1     0    0
2     1    0
3     2    0
4     3    0
5     4    0
6     0    1
7     1    1
8     2    1
9     3    1
10    4    1
11    0    2
12    1    2
13    2    2
14    3    2
15    4    2
16    0    3
17    1    3
18    2    3
19    3    3
20    4    3
21    0    4
22    1    4
23    2    4
24    3    4
25    4    4

直積の状態空間Sと状態ベクトルXを与えると、一番左の数字Yを返す関数を作った。

mmdp_X.to.Y = function(X,S) {
  mdp_check.state = function(S) {
    if( min(S)!=0 ) {
      stop("state space should begin with 0.")
    } else if( max(S)!=(length(S)-1) ) {
      stop("state space should end with N.")
    } else {
      return(TRUE)
    }
  }
  m = length(X)
  n = NULL
  for( i in c(1:m) ) {
    mdp_check.state(S[[i]])
    if( sum(X[i] %in% S[[i]]) == 0 ) {
      stop(paste("invalid first state. first state should be in {",paste(S.1, sep=" ", collapse=" "),"}."));
    }
    n = append(n,length(S[[i]]))
  }
  n = cumprod(n)
#  print(n)
#  print(X)
  y = n[-m]*X[-1]
  y = sum(y)+X[1]
  return(y+1)
}

これを使うには、直積を行う前の状態空間と、それぞれの次元の状態を次のように与える。

> mmdp_X.to.Y(c(4,3),list(c(0:4),c(0:4))) # 実行例
[1] 20

ここでリストオブジェクトを生成するために関数listを用いた。 状態はc(１次元目の状態, ２次元目の状態, …)のように与え、状態空間はlist(１次元目の状態空間, ２次元目の状態空間, …)のように与える。

またその直積空間の要素数をすべて数える関数も作った。

mmdp_X.to.Y.n = function(S) {
  mdp_check.state = function(S) {
    if( min(S)!=0 ) {
      stop("state space should begin with 0.")
    } else if( max(S)!=(length(S)-1) ) {
      stop("state space should end with N.")
    } else {
      return(TRUE)
    }
  }
  m = length(S)
  x.max = NULL
  for( i in c(1:m) ) {
    mdp_check.state(S[[i]])
    x.max = append(x.max,max(S[[i]]))
  }
  return(mmdp_X.to.Y(x.max,S))
}

この関数には状態空間のみを与える。

> mmdp_X.to.Y.n(list(c(0:4),c(0:4))) # 実行例
[1] 25

以上のような、直積空間の１次元への変換という考え方を用いれば、任意のコンポーネント数や任意の属性数のシステムのマルコフモデリングや、マルコフ決定過程の適用が可能になる。

これは各コンポーネントの劣化データから、markovchainパッケージの中の関数を用いて推定すれば良いため、省略する。

これは、すべての状態において、状態を新品0に戻す行動である。 そのような状態遷移行列は例えば、次の関数を用いて生成できる。

mmdp_create.replacement.matrix = function(S) {
  R.S = max(S)-min(S)
  n.S = length(S)
  if( R.S != n.S-1 ) {
    stop("state space is not regular and/or does not begin with 0.")
  }
  P = matrix(c(1,rep(0,n.S-1)),nrow=n.S,ncol=n.S,byrow=TRUE)
  rownames(P)=S
  colnames(P)=S
  return(P)
}

たとえば状態空間が{0,1,2,3,4}のシステムでの取替という行動は、次のように表される。

> mmdp_create.replacement.matrix(c(0:4))
  0 1 2 3 4
0 1 0 0 0 0
1 1 0 0 0 0
2 1 0 0 0 0
3 1 0 0 0 0
4 1 0 0 0 0

この表示された行列の意味を少し確認すること。

１次元ごとの状態遷移行列があるときに、直積空間上の状態遷移行列を作るには、それらを組み合わせる必要がある。その組み合わせを計算するために、次の関数を用意した。

mmdp_expand.P.2 = function(P) {
  n.1 = dim(P[[1]])[1]
  n.2 = dim(P[[2]])[1]
  P.Y = matrix(0,
               nrow=mmdp_X.to.Y.n(list(c(0:(n.1-1)),c(0:(n.2-1)))),
               ncol=mmdp_X.to.Y.n(list(c(0:(n.1-1)),c(0:(n.2-1)))))
  for( i in c(0:(n.1-1)) ) {
    for( j in c(0:(n.2-1)) ) {
      for( k in c(0:(n.1-1)) ) {
        for( l in c(0:(n.2-1)) ) {
          P.Y[mmdp_X.to.Y(c(i,j),list(c(0:(n.1-1)),c(0:(n.2-1)))),
              mmdp_X.to.Y(c(k,l),list(c(0:(n.1-1)),c(0:(n.2-1))))] = P[[1]][i+1,k+1]*P[[2]][j+1,l+1]
        }
      }
    }
  }
  return(P.Y)
}

この関数は、各次元の状態遷移が互いに独立であることを仮定して、同時確率を周辺確率の積で求める。

例題として劣化の状態遷移行列を用意する。

P.Dgr = matrix(c(
  9/10, 1/10, 0, 0, 0,
  0, 9/10, 1/10, 0, 0,
  0, 0, 9/10, 1/10, 0,
  0, 0, 0, 9/10, 1/10,
  0, 0, 0, 0, 1), nrow=5, ncol=5, byrow=TRUE)
rownames(P.Dgr) = c(0:4)
colnames(P.Dgr) = c(0:4)

この遷移行列に従うコンポーネントが二つあり、互いに独立に遷移する場合、２コンポーネントの同時の状態空間は、次のようになる。

expand.grid(c(0:4),c(0:4))

この空間での状態遷移行列の生成は、先ほど用意した関数mmdp_expand.P.2を用いて、次のように生成する。

> mmdp_expand.P.2(list(P.Dgr,P.Dgr)) # 実行例
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
 [1,] 0.81 0.09 0.00 0.00 0.00 0.09 0.01 0.00 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
 [2,] 0.00 0.81 0.09 0.00 0.00 0.00 0.09 0.01 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
 [3,] 0.00 0.00 0.81 0.09 0.00 0.00 0.00 0.09 0.01  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
 [4,] 0.00 0.00 0.00 0.81 0.09 0.00 0.00 0.00 0.09  0.01  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
 [5,] 0.00 0.00 0.00 0.00 0.90 0.00 0.00 0.00 0.00  0.10  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
 [6,] 0.00 0.00 0.00 0.00 0.00 0.81 0.09 0.00 0.00  0.00  0.09  0.01  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
 [7,] 0.00 0.00 0.00 0.00 0.00 0.00 0.81 0.09 0.00  0.00  0.00  0.09  0.01  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
 [8,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.81 0.09  0.00  0.00  0.00  0.09  0.01  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
 [9,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.81  0.09  0.00  0.00  0.00  0.09  0.01  0.00  0.00  0.00  0.00  0.00  0.00  0.00
[10,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.90  0.00  0.00  0.00  0.00  0.10  0.00  0.00  0.00  0.00  0.00  0.00  0.00
[11,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.81  0.09  0.00  0.00  0.00  0.09  0.01  0.00  0.00  0.00  0.00  0.00
[12,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.00  0.81  0.09  0.00  0.00  0.00  0.09  0.01  0.00  0.00  0.00  0.00
[13,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.00  0.00  0.81  0.09  0.00  0.00  0.00  0.09  0.01  0.00  0.00  0.00
[14,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.00  0.00  0.00  0.81  0.09  0.00  0.00  0.00  0.09  0.01  0.00  0.00
[15,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.00  0.00  0.00  0.00  0.90  0.00  0.00  0.00  0.00  0.10  0.00  0.00
[16,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.81  0.09  0.00  0.00  0.00  0.09  0.01
[17,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.81  0.09  0.00  0.00  0.00  0.09
[18,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.81  0.09  0.00  0.00  0.00
[19,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.81  0.09  0.00  0.00
[20,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.90  0.00  0.00
[21,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.90  0.10
[22,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.90
[23,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
[24,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
[25,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
      [,23] [,24] [,25]
 [1,]  0.00  0.00  0.00
 [2,]  0.00  0.00  0.00
 [3,]  0.00  0.00  0.00
 [4,]  0.00  0.00  0.00
 [5,]  0.00  0.00  0.00
 [6,]  0.00  0.00  0.00
 [7,]  0.00  0.00  0.00
 [8,]  0.00  0.00  0.00
 [9,]  0.00  0.00  0.00
[10,]  0.00  0.00  0.00
[11,]  0.00  0.00  0.00
[12,]  0.00  0.00  0.00
[13,]  0.00  0.00  0.00
[14,]  0.00  0.00  0.00
[15,]  0.00  0.00  0.00
[16,]  0.00  0.00  0.00
[17,]  0.01  0.00  0.00
[18,]  0.09  0.01  0.00
[19,]  0.00  0.09  0.01
[20,]  0.00  0.00  0.10
[21,]  0.00  0.00  0.00
[22,]  0.10  0.00  0.00
[23,]  0.90  0.10  0.00
[24,]  0.00  0.90  0.10
[25,]  0.00  0.00  1.00

同様に、片方は劣化、もう一方は交換するような行動の状態遷移行列は

mmdp_expand.P.2(list(P.Dgr,mmdp_create.replacement.matrix(c(0:4))))

あるいは２つ目の方を取り替えるのであれば

mmdp_expand.P.2(list(mmdp_create.replacement.matrix(c(0:4)),P.Dgr))

とする。

両方とも取替するのであれば

mmdp_expand.P.2(list(mmdp_create.replacement.matrix(c(0:4)),mmdp_create.replacement.matrix(c(0:4))))

で生成される。

２つのマルコフ決定過程それぞれの報酬あるいは費用の組み合わせも、１次元に射影した直積空間に変換できる。ここでは、２つの報酬あるいは費用は、それぞれの報酬あるいは費用の和と仮定している。

mmdp_expand.R.2 = function(R) {
  n.1 = length(R[[1]])
  n.2 = length(R[[2]])
  R.Y = rep(0,mmdp_X.to.Y.n(list(c(0:(n.1-1)),c(0:(n.2-1)))))
  for( i in c(0:(n.1-1)) ) {
    for( j in c(0:(n.2-1)) ) {
      R.Y[mmdp_X.to.Y(c(i,j),list(c(0:(n.1-1)),c(0:(n.2-1))))] = R[[1]][i+1]+R[[2]][j+1]
    }
  }
  return(R.Y)
}

この関数も次のようにリストを引数に持たせて使う。

> mmdp_expand.R.2(list(c(1,2,3,4,5),c(1,2,3,4,5))) # 実行例
 [1]  2  3  4  5  6  3  4  5  6  7  4  5  6  7  8  5  6  7  8  9  6  7  8  9 10

以上の関数を用いると、２次元のマルコフ決定過程を１次元に射影できる。

前週の実験で推定した製品を２つ使うシステムを考えると、直積空間での劣化の状態遷移行列は次のようになる。

P.Dgr.Dgr = mmdp_expand.P.2(list(P.Dgr,P.Dgr))

ひとつ目を交換し、２つ目はそのまま劣化させる行動は、次のような状態遷移行列を持つ。

P.Rpl.Dgr = mmdp_expand.P.2(list(mmdp_create.replacement.matrix(c(0:(dim(P.Dgr)[1]-1))),
                                    P.Dgr))

逆に、ひとつ目はそのまま劣化させ、２つ目を交換する行動は、次のような状態遷移行列を持つ。

P.Dgr.Rpl = mmdp_expand.P.2(list(P.Dgr,
                                    mmdp_create.replacement.matrix(c(0:(dim(P.Dgr)[1]-1)))))

両方とも交換する行動の状態遷移行列も作れる。

P.Rpl.Rpl = mmdp_expand.P.2(list(
                                    mmdp_create.replacement.matrix(c(0:(dim(P.Dgr)[1]-1))),
                                    mmdp_create.replacement.matrix(c(0:(dim(P.Dgr)[1]-1)))))

最後にこれらを、一つの配列に納めて、MDPtoolboxで使う準備をする。

P = array(0, dim=c(dim(P.Dgr.Dgr)[1],dim(P.Dgr.Dgr)[2],4))
P[,,1] = P.Dgr.Dgr
P[,,2] = P.Rpl.Dgr
P[,,3] = P.Dgr.Rpl
P[,,4] = P.Rpl.Rpl

費用も同様に、コンポーネントごとに用意して、それらを一つの配列に収める。

C.Opr = c(0,0,0,0,2000)
C.Rpl = c(150,150,150,150,150)
C.Opr.Opr = mmdp_expand.R.2(list(C.Opr,C.Opr))
C.Rpl.Opr = mmdp_expand.R.2(list(C.Rpl,C.Opr))
C.Opr.Rpl = mmdp_expand.R.2(list(C.Opr,C.Rpl))
C.Rpl.Rpl = mmdp_expand.R.2(list(C.Rpl,C.Rpl))
Cost = cbind(C.Opr.Opr,C.Rpl.Opr,C.Opr.Rpl,C.Rpl.Rpl)
R = -Cost

組み合わせのコストが足し算と異なる場合には、各自で設定する必要がある。

あとは実行するだけ。 例によって準備をする。

install.packages("MDPtoolbox")
library(MDPtoolbox)

そして、上の例の最適方策を価値反復アルゴリズムによって求めると、次のような結果を得る。

> mdp_value_iteration(P,R,0.95)
[1] "MDP Toolbox: iterations stopped, epsilon-optimal policy found"
$V
 [1]  -64.03890  -81.67123 -108.58531 -149.66496 -212.36385  -81.67123  -99.30356 -126.21764 -167.29729 -229.99618 -108.58531
[12] -126.21764 -153.13171 -194.21137 -256.91025 -149.66496 -167.29729 -194.21137 -235.29102 -297.98991 -212.36385 -229.99618
[23] -256.91025 -297.98991 -360.68879

$policy
 [1] 1 1 1 1 2 1 1 1 1 2 1 1 1 1 2 1 1 1 1 2 3 3 3 3 4

$iter
[1] 76

$time
Time difference of 0.04542994 secs

$epsilon
[1] 0.01

$discount
[1] 0.95

生成された最適方策policyがどのような意味を持つかは、次のように状態空間と照らしてみる必要がある。

> optimal.policy = mdp_value_iteration(P,R,0.95) # 実行例
> cbind(expand.grid(c(0:4),c(0:4)),optimal.policy$policy,optimal.policy$V) # 実行例
   Var1 Var2 optimal.policy$policy optimal.policy$V
1     0    0                     1        -64.03890
2     1    0                     1        -81.67123
3     2    0                     1       -108.58531
4     3    0                     1       -149.66496
5     4    0                     2       -212.36385
6     0    1                     1        -81.67123
7     1    1                     1        -99.30356
8     2    1                     1       -126.21764
9     3    1                     1       -167.29729
10    4    1                     2       -229.99618
11    0    2                     1       -108.58531
12    1    2                     1       -126.21764
13    2    2                     1       -153.13171
14    3    2                     1       -194.21137
15    4    2                     2       -256.91025
16    0    3                     1       -149.66496
17    1    3                     1       -167.29729
18    2    3                     1       -194.21137
19    3    3                     1       -235.29102
20    4    3                     2       -297.98991
21    0    4                     3       -212.36385
22    1    4                     3       -229.99618
23    2    4                     3       -256.91025
24    3    4                     3       -297.98991
25    4    4                     4       -360.68879

この方策は「状態が４になったのみ方を取替し、両方とも４になったら両方とも取替する」という方策であると、読み取れるだろうか。