推移行列・行列の指数関数
- Rmd
推移行列・行列の指数関数: Transition Matrix and Matrix Exponential
- 作者: ryamada
- 発売日: 2016/07/21
- メディア: Kindle版
- この商品を含むブログ (1件) を見る
--- title: "Transition Matrix 状態推移行列" author: "ryamada" date: "Friday, July 22, 2016" output: html_document --- # Discrete time formula 離散時間の式 $$ x(t+1) = -y(t)\\ y(t+1) = x(t) $$ Calculate values step-by-step. ```{r} n.t <- 20 # No. steps # Initial values x0 <- 2 y0 <- 1.5 xy.init <- c(x0,y0) xy <- matrix(0,n.t,2) xy[1,] <- xy.init for(i in 2:n.t){ xy[i,1] <- -xy[i-1,2] xy[i,2] <- xy[i-1,1] } par(mfcol=c(1,2)) matplot(xy,type="l") plot(xy,type="l") par(mfcol=c(1,1)) ``` # Matrix 行列 $$ x(t+1) = 0 \times x(t) + (-1) \times y(t)\\ y(t+1) = 1 \times x(t) + 0 \times y(t) $$ Calculate with the matrix. ```{r} M <- matrix(c(0,-1,1,0),2,2,byrow=TRUE) M xy.2 <- xy for(i in 2:n.t){ xy.2[i,] <- M %*% xy.2[i-1,] } par(mfcol=c(1,2)) matplot(xy.2,type="l") plot(xy.2,type="l") par(mfcol=c(1,1)) ``` # Periodic movement -> Rotation 周期的運動なので回転と見做してみる $$ x(t+1) = \cos{\theta} \times x(t) + (-\sin{\theta}) \times y(t)\\ y(t+1) = \sin{\theta} \times x(t) + \cos{\theta} \times y(t) $$ ```{r} theta <- pi/11 M <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),2,2,byrow=TRUE) M ``` # Eigen value decomposition and power of matrix 固有値分解と行列の冪 ```{r} eigen.out <- eigen(M) eigen.out[[1]] eigen.out[[2]] ``` $$ M = V S U\\ VU = UV =I\\ S = diag(\lambda) $$ ```{r} V <- eigen.out[[2]] S <- diag(eigen.out[[1]]) V U <- solve(V) V %*% U U %*% V S V %*% S %*% U V %*% S %*% U - M ``` $$ S^k = diag(\lambda^k) $$ ```{r} S %*% S # S^2 diag(eigen.out[[1]]^2) M %*% M (V %*% S %*% U) %*% (V %*% S %*% U) # V %*% S %*% (U %*% V) %*% S %*% U # V %*% S %*% I %*% S %*% U # V %*% S %*% S %*% U # V %*% S^2 %*% U V %*% (diag(eigen.out[[1]]^2)) %*% U ``` $$ M^k = V S^k U $$ ```{r} # Power of matrix my.pow.mat <- function(M,k){ eigen.out <- eigen(M) V <- eigen.out[[2]] U <- solve(V) lambdas <- eigen.out[[1]] + 0 * 1i # Force to handle the values as complex ret <- V %*% diag(lambdas^k) %*% U return(ret) } my.pow.mat(M,1) my.pow.mat(M,2) M %*% M my.pow.mat(M,0) solve(M) my.pow.mat(M,-1) ``` Imaginary components are 0, so return real components only. 虚部は0になるので、返り値を実部のみにする。 ```{r} my.pow.mat <- function(M,k){ eigen.out <- eigen(M) V <- eigen.out[[2]] U <- solve(V) lambdas <- eigen.out[[1]] + 0 * 1i ret <- V %*% diag(lambdas^k) %*% U return(Re(ret)) } my.pow.mat(M,1) my.pow.mat(M,2) M %*% M my.pow.mat(M,0) solve(M) my.pow.mat(M,-1) ``` k can be real values. kは実数でもよい。 Calculate valuew using power of the matrix. 行列の冪を使って計算する ```{r} n.t <- 1000 ks <- seq(from=0,to=50,length=n.t) xy.3 <- matrix(0,n.t,2) for(i in 1:n.t){ xy.3[i,] <- my.pow.mat(M,ks[i]) %*% xy.init # Initial values } par(mfcol=c(1,2)) matplot(xy.3,type="l") plot(xy.3,type="l") par(mfcol=c(1,1)) ``` ```{r} n.t <- 1000 ks <- seq(from=0,to=50,length=n.t) n.v <- 2 # No. variables might be changed # Random transition matrix M <- matrix(rnorm(n.v^2),n.v,n.v) M xy.4 <- matrix(0,n.t,n.v) xy.init <- rnorm(n.v) for(i in 1:n.t){ xy.4[i,] <- my.pow.mat(M,ks[i]) %*% xy.init } par(mfcol=c(1,2)) matplot(xy.4,type="l") plot(xy.4,type="l") par(mfcol=c(1,1)) ``` Various patterns for n.v =2. 色々な変化(n.v=2)。 ```{r} n.iter <- 10 for(ii in 1:n.iter){ n.t <- 1000 ks <- seq(from=0,to=50,length=n.t) n.v <- 2 M <- matrix(rnorm(n.v^2),n.v,n.v) M xy.4 <- matrix(0,n.t,n.v) xy.init <- rnorm(n.v) for(i in 1:n.t){ xy.4[i,] <- my.pow.mat(M,ks[i]) %*% xy.init } par(mfcol=c(1,2)) matplot(xy.4,type="l") plot(xy.4,type="l") par(mfcol=c(1,1)) } ``` Change n.v to 3. 変数の数を3に変える。 ```{r} n.iter <- 10 for(ii in 1:n.iter){ n.t <- 1000 ks <- seq(from=0,to=20,length=n.t) n.v <- 3 M <- matrix(rnorm(n.v^2),n.v,n.v) M xy.4 <- matrix(0,n.t,n.v) xy.init <- rnorm(n.v) for(i in 1:n.t){ xy.4[i,] <- my.pow.mat(M,ks[i]) %*% xy.init } par(mfcol=c(1,2)) matplot(xy.4,type="l") plot(as.data.frame(xy.4),cex=0.1,pch=20) par(mfcol=c(1,1)) } ``` # System of Differential Equation 連立微分方程式 $$ \begin{pmatrix} \frac{dx_1}{dt} \\ \frac{dx_2}{dt} \end{pmatrix} = \begin{pmatrix} a & b \\ c & d \end{pmatrix} \begin{pmatrix} x_1\\ x_2 \end{pmatrix} $$ The solution of this system is known as, この連立微分方程式の解は以下であると知られている。 ! Study "differential equations and matrix" yourself. $$ \begin{pmatrix} x_1(t)\\ x_2(t) \end{pmatrix} = e^{A t} \begin{pmatrix} x_1(0)\\ x_2(0) \end{pmatrix}\\ A = V diag(\lambda_A) U\\ \Sigma(t) = diag(e^{\lambda_A t})\\ e^{A t} = V \Sigma(t) U $$ This is similar to the previous sections' expressions. この式は、前項で扱ってきた式によく似ている。 $$ \begin{pmatrix} x_1(t+1)\\ x_2(t+1) \end{pmatrix} = M \begin{pmatrix} x_1(t)\\ x_2(t) \end{pmatrix}\\ \begin{pmatrix} x_1(t)\\ x_2(t) \end{pmatrix} = M(t) \begin{pmatrix} x_1(0)\\ x_2(0) \end{pmatrix}\\ M = V diag(\lambda_M) U\\ M(t) = V \Sigma(t) U\\ \Sigma(t) = diag(\lambda_M^t)\\ \lambda_M = e^{\lambda_A} $$ ```{r} n.v <- 2 theta <- pi/11 M <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),byrow=TRUE,2,2) M[1,1] <- M[1,1] + 0 M ``` Transform into differential eqations. ```{r} my.mat.exp <- function(M,t){ eigen.out <- eigen(M) V <- eigen.out[[2]] U <- solve(V) lambda <- eigen.out[[1]] + 0 * 1i ret <- V %*% diag(exp(lambda*t)) %*% U return(Re(ret)) } my.mat.pow2exp <- function(M){ eigen.out <- eigen(M) V <- eigen.out[[2]] U <- solve(V) lambda <- eigen.out[[1]] + 0 * 1i ret <- V %*% diag(log(lambda)) %*% U return(Re(ret)) } my.mat.exp2pow <- function(M){ return(my.mat.exp(M,1)) } A <- my.mat.pow2exp(M) ``` Transision matrix for a unit time interval: M Differential equation matrix for infinitesimal time interval : A 単位時間ごとの状態推移行列:M 極小時間ごとの微分方程式行列:A ```{r} M A ``` ```{r} n.t <- 30 xy.5 <- xy.6 <- matrix(0,n.t,n.v) xy.init <- rnorm(n.v) xy.5[1,] <- xy.init xy.6[1,] <- my.mat.exp(A,ks[1]) %*% xy.init for(i in 2:n.t){ xy.5[i,] <- M %*% xy.5[i-1,] xy.6[i,] <- my.mat.exp(A,i) %*% xy.init } par(mfcol=c(1,2)) matplot(xy.5,type="l") matplot(xy.6,type="l") par(mfcol=c(1,1)) ``` # Transision Probability Matrix 推移確率行列 Multiple n states. 複数(n個)の状態 Probability vector is a vector with n elements; 確率ベクトルはn個の要素を持ち、 $$ P(t) = (p_1,...,p_n)\\ p_i >= 0\\ \sum_{i=1}^n p_i= 1 $$ ## Hospital 病院 Assume four states: (1) hospitalization for tests (T), (2) sick (S), (3) Very sick (V), and (4) vacant/blank (B). B stands for vacancy due to discharge including death-discharge. 病院のベッドの状態が4つあるとする。検査入院(T)、軽症入院(S)、重症入院(S)、空ベッド(B)。 After a unit time period (a day), patients in state T change their status to (T,S,V,B) with probability (0.5,0.04,0.01,0.45); average hospital stays for tests are approximately 2 days and very likely to be discharged without clinical events. 単位時間(1日)ごとに、各状態の患者はその状態を維持するか変えるかする。 例えば、検査入院の人は、次の日には、(0.5,0.04,0.01,0.45)の確率で(T,S,V,B)の状態に移ると言う。半数はそのまま検査入院、4%は軽症状態になり、1%は重症状態になる。残りは退院する(ふつうは検査が終わって退院する)。 This is written as T : (0.5,0.04,0.01,0.45). それを、T : (0.5,0.04,0.01,0.45)と書くことにすれば、 In the same way, S : (0,0.7,0.1,0.2) V : (0,0.05,0.8,0.15) B : (0.6,0.3,0.1,0) ```{r} t <- c(0.5,0.04,0.01,0.45) s <- c(0,0.7,0.1,0.2) v <- c(0,0.05,0.8,0.15) b <- c(0.6,0.3,0.1,0) M <- cbind(t,s,v,b) apply(M,2,sum) # should be 1 n.t <- 30 x.init <- c(0,0,0,1) x <- matrix(0,n.t,4) x[1,] <- x.init for(i in 2:n.t){ x[i,] <- M %*% x[i-1,] # x[i,] <- my.mat.pow(M,i) %*% x.init } matplot(x,type="l") apply(x,1,sum) ``` How do you reduce the fraction of B. どうやって空きベッド率を減らすのがよいか？
