# 推移行列・行列の指数関数

• Rmd

---
title: "Transition Matrix 状態推移行列"
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.

{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.

{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.

{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

{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;

$$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.

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.

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.

どうやって空きベッド率を減らすのがよいか？