ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2016-07-26 Rで探し物

[]Rで探し物

2016-07-25 R TaskViewを使って独り立ち

[][][]R TaskViewを使って独り立ちする

  • ある程度Rが使えるようになったら、あとは自分がやりたいことをするための調べ物ができるようになることが大事
  • 短期集中セミナーである程度Rが使えるようになったら、調べ物のやり方について知ることが大事
  • そんな短期集中セミナーの最終回のテキスト(Rmd フォーマット)
---
title: "TimeSeries時系列"
author: "ryamada"
date: "2016年7月24日"
output: html_document
---

# R for Time Series Analyses 時系列解析のためのR

When you are new to a type of analysis, "time series analysis" for example, you may want to check packages useful for it.

新しいタイプの解析、たとえば時系列解析、をしようとするときに、有用なパッケージを探したいと思うかもしれない。

Task view is one place to look for. "Task view"に探しに行くこともできる。

Goocle with "CRAN time series task view".
"CRAN time series task view"でグーグル。

https://cran.r-project.org/web/views/TimeSeries.html


# List topics handled in the article 記事で扱っているトピックスを確認

- Basics
- Times and Dates
- Time Series Classes
- Forecasting and Univariate Modeling
- Frequency analysis
- Decomposition and Filtering
- Seasonality
- Stationarity, Unit Roots, and Cointegration
- Nonlinear Time Series Analysis
- Dynamic Regression Models
- Multivariate Time Series Models
- Continuous time models
- Resampling
- Time Series Data
- Miscellaneous

# Understand Time Series Analyses Based on the List of Topics トピックスリストに基づいて時系列解析を理解する

## Data handling
* Basics
    + Regularly spaced time series (using numeric time stamps)
    + 規則的な時間間隔系列。数値でタイムスタンプ
    + IRRegularly spaced data needs something special
    + 不規則な時間間隔データには特別な何かが必要
* Time Series Classes
    + ts class is the base and there are many other classes. Grab the ts class first and then select one of others if necessary
    + tsクラスを基本にたくさんのクラスがあるから、tsの概要を理解した上で必要に応じて選ぶ

## Our task is what kind of time-series analysis and how? 自分のやりたいことはどういうタイプの時系列解析で、どうやるか?
* Tasks handled 扱う仕事
    + Forecasting and Univariate Modeling 予測・1変数
    + Frequency analysis 周波数解析
    + Decomposition and Filtering 分解・フィルタリング
    + Seasonality 季節性変動
    + Stationarity, Unit Roots, and Cointegration 定常性…?
* How to analyze どう解析するか
    + Nonlinear Time Series Analysis
    + Linear analysis is the basic and NONlinear is special
    + Dynamic Regression Models
    + Regular regression analysis is "Non-dynamic"
    + Multivariate Time Series Models
    + Mono-variate is regular and MULTI-variate is special
    + Continuous time models
    + NON-continuous = DISCRETE is regular and Continuous is special
    


# Learn with examples Examples を使って学ぶ

## ts

```{r}
help(ts)
```

```{r}
#example(ts)
```

```{r}
require(graphics)

ts(1:10, frequency = 4, start = c(1959, 2)) # 2nd Quarter of 1959
print( ts(1:10, frequency = 7, start = c(12, 2)), calendar = TRUE)
# print.ts(.)
## Using July 1954 as start date:
gnp <- ts(cumsum(1 + round(rnorm(100), 2)),
          start = c(1954, 7), frequency = 12)
plot(gnp) # using 'plot.ts' for time-series plot

## Multivariate
z <- ts(matrix(rnorm(300), 100, 3), start = c(1961, 1), frequency = 12)
class(z)
head(z) # as "matrix"
plot(z)
plot(z, plot.type = "single", lty = 1:3)

## A phase plot:
plot(nhtemp, lag(nhtemp, 1), cex = .8, col = "blue",
     main = "Lag plot of New Haven temperatures")

```

Look inside of objects. オブジェクトを見てみる
```{r}
gnp
head(gnp)
str(gnp) # structure
class(gnp)
```

## Decomposition and Filtering

```
Filters and smoothing : filter() in stats provides autoregressive and moving average linear filtering of multiple univariate time series. The robfilter package provides several robust time series filters, while mFilter includes miscellaneous time series filters useful for smoothing and extracting trend and cyclical components. smooth() from the stats package computes Tukey's running median smoothers, 3RS3R, 3RSS, 3R, etc. sleekts computes the 4253H twice smoothing method.
```

```{r}
help(filter)
# example(filter)
x <- rnorm(100)
x <- cumsum(x)
par(mfrow=c(2,3))
plot(x,type="l",main="original")
tmp.x <- x
for(i in 2:6){
  tmp.x <- filter(tmp.x,rep(1,3)/3) 
  plot(tmp.x,col=i,type="l",main=paste("",i-1))
}
```

```{r}
help(decompose)
# example(decompose)
require(graphics)
plot(co2)
m <- decompose(co2)
m$figure
plot(m)

## example taken from Kendall/Stuart
x <- c(-50, 175, 149, 214, 247, 237, 225, 329, 729, 809,
       530, 489, 540, 457, 195, 176, 337, 239, 128, 102, 232, 429, 3,
       98, 43, -141, -77, -13, 125, 361, -45, 184)
x <- ts(x, start = c(1951, 1), end = c(1958, 4), frequency = 4)
m <- decompose(x)
## seasonal figure: 6.25, 8.62, -8.84, -6.03
round(decompose(x)$figure / 10, 2)
```
```{r}
x <- ts(cumsum(rnorm(1000)),frequency=12) # No seasonality
m <- decompose(x)
plot(m)
str(m)
```



# Optional

## Packages and Task Views パッケージとTask View

See https://mran.revolutionanalytics.com/rpackages/ 

You can download ALL PACKAGES listed in a task view at once, (if you really want to do so).

あるTask viewに扱われているパッケージを一括でダウンロードすることもできる(もし本当にそうしたいのなら)## Search tools 探し物

http://www.r-bloggers.com/21-r-navigation-tools/ 

2016-07-24 乱数とシミュレーション

[][][]乱数とシミュレーション

  • Rmd
---
title: "RandomValues and Simulation"
author: "ryamada"
date: "2016年7月23日"
output: html_document
---

# Generate Random Values 乱数を作る

## Uniform distribution 一様分布

More numbers, more uniform. 数が増えると一様に近づく。
Generate many numbers and use a part of them for small sample size exapmles. 一度にたくさんの乱数を作っておいて、その一部を少サンプルの例として使う。

```{r}
n <- 10^6
r <- runif(n)
par(mfrow=c(2,3))
for(i in 1:6){
  tmp <- r[1:(10^i)]
  hist(tmp)
}
par(mfrow=c(1,1))
```

Sample mean converges to 0.5. 標本平均は0.5に近づく。

```{r}
n <- 10^4
r <- runif(n)
means <- rep(0,n)
for(i in 1:n){
  tmp <- r[1:i]
  means[i] <- mean(tmp)
}
plot(means,type="l",ylim=c(0,1))
```

Some functions work on vectors quicker than loop.
ベクトルに適用する関数には、ループを回すより速いものがある。

```{r}
# quicker
cmsm <- cumsum(r)
means.2 <- cmsm/(1:n)
plot(means.2,type="l",ylim=c(0,1))
```

## Normal distribution 正規分布
```{r}
n <- 10^4
hist(rnorm(n))
```

## Multidimensional 多次元

## (Square, Cube,...) vs. (Circle, Sphere,...) (正方形、立方体、…) 対 (円、球、…)

Square uniform distribution in 2D 2次元で正方領域一様分布。

Two-dimensional normal distribution. 2次元正規分布

```{r}
d <- 2
n <- 10^4
xy1 <- matrix(runif(n*d),ncol=d)
plot(xy1)

xy2 <- matrix(rnorm(n*d),ncol=d)
plot(xy2)
```

3D cubic. 3次元立方体。3D normal. 3次元正規分布
。
```{r setup}
library(knitr)
library(rgl)
knit_hooks$set(rgl = hook_rgl)
d <- 3
xy1 <- matrix(runif(n*d),ncol=d)
```

```{r, rgl=TRUE}
plot3d(xy1)
```
```{r, rgl=TRUE}
xy2 <- matrix(rnorm(n*d),ncol=d)
plot3d(xy2)
```

## Directionally Uniform 方向一様

Directional evaluation needs directinally uniform distribution.
方向についての評価では、方向に関する一様分布が必要になる。

Normal distribution is equal to all directions in any dimension.

正規分布は次元によらず、全方向について平等。

For 2D, uniform distribution of angle can be used.

2次元では、角度に関する一様分布も使える

```{r}
d <- 2
xy2 <- matrix(rnorm(n*d),ncol=d)
L <- sqrt(apply(xy2^2,1,sum))
xy2. <- xy2/L
plot(xy2.)
theta <- runif(n) * 2 * pi
xy3. <- cbind(cos(theta),sin(theta))
plot(xy3.)
```

3D. 3次元。

```{r, rgl=TRUE}
d <- 3
xy2 <- matrix(rnorm(n*d),ncol=d)
L <- sqrt(apply(xy2^2,1,sum))
xy2. <- xy2/L
plot3d(xy2.)
# ...
```

1D directionally uniform distribution...
1次元の方向一様分布とは?

```{r}
d <- 1
xy2 <- matrix(rnorm(n*d),ncol=d)
L <- sqrt(apply(xy2^2,1,sum))
xy2. <- xy2/L
plot(xy2.)
xy3. <- sample(c(-1,1),n,replace=TRUE)
plot(xy3.)
```

## Disk (circle inside) 円板

Disk is a peripheral circle with inside.
円板は円周と内部。

Uniform distribution of disk. 円板の一様分布。

Take random values only satisfying the rule. 条件に合う乱数のみを採用する。

```{r}
d <- 2
xy1 <- matrix(runif(n*d)*2-1,ncol=d)
L <- sqrt(apply(xy1^2,1,sum))
disk <- which(L<=1)
xy1. <- xy1[disk,]
plot(xy1.)
length(disk) / n
```

3D 3次元。

Fraction meeting the rule get lower.
採用割合が下がる。
```{r, rgl=TRUE}
d <- 3
xy1 <- matrix(runif(n*d)*2-1,ncol=d)
L <- sqrt(apply(xy1^2,1,sum))
disk <- which(L<=1)
xy1. <- xy1[disk,]
plot3d(xy1.)
length(disk) / n
```

10D's fraction is really low. 10次元の合格乱点割合は本当に低い。
```{r, rgl=TRUE}
d <- 10
xy1 <- matrix(runif(n*d)*2-1,ncol=d)
L <- sqrt(apply(xy1^2,1,sum))
disk <- which(L<=1)
xy1. <- xy1[disk,]
plot3d(xy1.)
length(disk)
length(disk) / n
```

## Discrete sampling 離散サンプリング

{0,1} even. {0,1}均等。
```{r}
n <- 10^4
s1 <- sample(0:1,n,replace=TRUE)
table(s1)
```

{0,1} with specific fraction. 割合を指定して{0,1}。

```{r}
n <- 10^4
p <- c(0.3,0.7)
s1 <- sample(0:1,n,replace=TRUE,prob=p)
table(s1)
```

n categories. カテゴリ数n。

```{r}
n <- 10^4
p <- c(0.5,0.2,0.1,0.1,0.1)
s1 <- sample(1:5,n,replace=TRUE,prob=p)
table(s1)
```

Number of samples per category can be generated with multinominal random generator. カテゴリ別のサンプル数のみを発生するなら多項分布乱数発生関数でできる。

```{r}
rmultinom(1,n,prob=p)
rmultinom(3,n,prob=p)
```

## Shuffle シャッフルする

sample() function can shuffle all, take out a part with or without replacement.

sample()関数は全要素の入れ替え、一部要素の取り出し(重複サンプリングについて、ありの場合となしの場合と)ができる。
```{r}
h <- c(1,4,5,7,9)
sample(h)
sample(h,replace=TRUE)
sample(h,3)
sample(h,10,replace=TRUE)
sample(h,4,replace=TRUE)
```

## Seeds シード

Seeds control random value generators. シードによって乱数発生関数をコントロールする。

```{r}
runif(3)
runif(3)
set.seed(12345)
runif(3)
runif(3)
set.seed(12345)
runif(3)
runif(3)
set.seed(12345)
runif(3)
set.seed(12345)
runif(3)
runif(3)
set.seed(12345)
runif(6)

set.seed(12346)
runif(3)
```
## Random generating functions 乱数発生関数のいろいろ

```{r}
help(Distributions)
```
dxxxx: probability density function 確率密度
rxxxx: random value generator 乱数発生

# Simulation シミュレーション

## Damages pile up 傷害が積み重なる

Assume a cell will get damages periodically with probability p per a day.
When the number of damaging events reach q, the cell dies.
Generate distribution of life spans simulationally.

細胞はダメージを受けるが、そのダメージの生起確率は1日当たりpだと云う。ダメージイベントの回数がqになると細胞は死ぬという。
細胞の生存期間の分布をシミュレーションで作れ。

```{r}
n <- 10^4
p <- 0.05
q <- 20
days <- 10^3
x <- matrix(sample(0:1,n*days,replace=TRUE,prob=c(1-p,p)),ncol=n)
#image(x)
x.cumsum <- apply(x,2,cumsum)
matplot(x.cumsum[,1:100],type="l")
abline(h=q)
```
```{r}
plot(x.cumsum[,1])
abline(h=q)
which(x.cumsum[,1]==q)
min(which(x.cumsum[,1]==1))
res <- apply(x.cumsum,2,function(x){min(which(x==q))})
hist(res)
plot(sort(res))
```

## Differential equations with random errors ランダムエラー付きの微分方程式

$$
\frac{dx}{dt} = p x - q xy\\
\frac{dy}{dt} = r xy - s y
$$

```{r}
n.t <- 10000
p <-3
q <- 2
r <- 4
s <- 1

x.init <- c(1,1)
dt <- 10^(-3)

x <- matrix(0,n.t,2)
x[1,] <- x.init
for(i in 2:n.t){
  dxdt <- dt*(p *x[i-1,1] - q * x[i-1,1]*x[i-1,2])
  dydt <- dt*(r * x[i-1,1] * x[i-1,2] - s * x[i-1,2])
  x[i,1] <- x[i-1,1] + dxdt
  x[i,2] <- x[i-1,2] + dydt
}
matplot(x,type="l")
plot(x,type="l")
```

```{r}
x <- matrix(0,n.t,2)
x[1,] <- x.init
for(i in 2:n.t){
  dxdt <- dt*(p *x[i-1,1] - q * x[i-1,1]*x[i-1,2])
  dydt <- dt*(r * x[i-1,1] * x[i-1,2] - s * x[i-1,2])
  x[i,1] <- x[i-1,1] + dxdt*(1+rnorm(1)*1)
  x[i,2] <- x[i-1,2] + dydt*(1+rnorm(1)*1)
}
matplot(x,type="l")
plot(x,type="l")
```

2016-07-23 多次元正規分布を作る

[][]多次元正規分布を作る

# 多次元正規分布乱数を作る関数を作っておく
library(MASS)
library(lqmm)
# 適当に対称行列を作って正定値行列を作る関数も作っておく
# 多次元正規分布は中心座標と分散共分散行列とで決まる
# 分散共分散行列のために正定値行列が必要
my.make.positive.definite <- function(n.v){# n.vは変数の個数
  sigma <- matrix(rnorm(n.v^2),n.v,n.v)
  sigma <- sigma + t(sigma)
  sigma <- make.positive.definite(sigma)
  return(sigma)
}
# n.sampleは発生乱点数、n.vは変数の個数(次元)
# mは中心座標、sigmaはn.v x n.v 正定値行列
my.mvrnorm <- function(n.sample,n.v,m=NULL,sigma=NULL){
  if(is.null(m)){
    m <- rnorm(n.v)
  }
  if(is.null(sigma)){
    sigma <- my.make.positive.definite(n.v)
  }
  return(mvrnorm(n.sample,m,sigma))
}


# サンプル数
n.s <- 1000

# 説明変数
## 説明変数の数
n.v <- 5
## 多次元正規分布とする
M <- rnorm(n.v) # 重心
S <- my.make.positive.definite(n.v) # Sigma行列(分散共分散行列のようなもの)
X <- my.mvrnorm(n.s,n.v,M,S)
plot(as.data.frame(X))

library(rgl)


tmp.X <- X[,1:3]
tmp.X <- rbind(tmp.X,rep(max(tmp.X),3),rep(min(tmp.X),3)) # 3軸スケールをそろえるため

plot3d(tmp.X)

# sigma行列を適当に与えて、多次元正規分布(ただしゆがんでいる)を作った

# 点の数をだんだん増やす
n.ss <- c(10,25,50,75,100,250,500,750,1000,5000,10000)

diff.M <- rep(0,length(n.ss))
for(i in 1:length(n.ss)){
	n.s <- n.ss[i]
	X <- my.mvrnorm(n.s,n.v,M,S)

	# この乱点のラグビーボール状態を調べるのは、PCA
	# 標本重心
	M.hat <- apply(X,2,mean)
	# 真値と比べる
	M.hat
	M

	pca.out <- prcomp(X)
	# 長軸からsdが並んで出る
	# sdを対角成分にした行列を作り

	K <- diag(pca.out$sdev)
	# 回転行列も出ているから、それで回す
	R <- pca.out$rotation

	P <- R %*% K

	Q <- P %*% t(P) # var-covとは各軸の内積のようなもの
	S

	diff.M[i] <- sum((Q-S)^2)

}

plot(log(n.ss),diff.M,type="l")

2016-07-22 推移行列・行列の指数関数

[][]推移行列・行列の指数関数

  • Rmd
---
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.

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