Untitled
状態推移 Status-transition
基本 Basics
- 確率的に何かが起きる Consider a stochastic process
- 次に何が起きるかは「今」の影響を受けるが、「前」の影響を受けない “Next event” depends on “Now” but not on “Past”.
簡単な例
ただの確率事象 A simple stochastic process
ある確率で、機能性細胞が死滅するとする。その確率はpである。初め、細胞は全部でn個あった。生存細胞数をシミュレーションする。
Assume that cells die stochasitically. The probability of death is p. Simulate for n cells at the begininng.
n <- 1000
p <- 0.1
s.or.d <- sample(0:1,n,replace=TRUE,prob=c(1-p,p))
s.or.d[1:100]## [1] 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 1
## [36] 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [71] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0sum(s.or.d) # 死細胞数 No. dead cells## [1] 105sum(s.or.d)/n # 死細胞割合 Fraction of dead cells## [1] 0.105増えたり減ったり Increase and decrease
細胞は確率的に死ぬが、生き残った細胞は増殖する。死亡確率をp、増殖(１細胞が２細胞になる)確率をqとする。初め、細胞は全部でn個あった。
Cells may die but living cells will divide into two. Both are stochastic, with probabilities p and q, respectively. At the begininng there are n cells.
いくつの細胞が死に、いくつの細胞が分裂するかは、その時々の生細胞数に依存するので、「今」の影響を受けている。 しかしながら、一段階前に何個の細胞が生きていたかには影響を受けていない。 No. death/division depend on the No. cells “NOW”, but not on “One period previous”.
p <- 0.05
q <- 0.03
weeks <- 1:300
# No. living cell from 1st through 100th week
n.sv.1 <- rep(0,length(weeks))
n.sv.1[1] <- 100
for(i in 2:length(weeks)){
# death
tmp1 <- sample(0:1,n.sv.1[i-1],replace=TRUE,prob=c(1-p,p))
n.death <- sum(tmp1)
# No. living cells for now
tmp2 <- n.sv.1[i-1] - n.death
# Cell division
tmp3 <- sample(0:1,tmp2,replace=TRUE,prob=c(1-q,q))
# No. new born cells
tmp4 <- sum(tmp3)
# No. living cells for now
n.sv.1[i] <- tmp2 + tmp4
}
plot(n.sv.1)
加速する Acceration
生きている細胞がm 個のとき、1/mの仕事をしなくてはならず、仕事量が多いと、それに応じて死亡率が上がるとする。分裂確率は仕事割合によって変わらないとする。
When there are m cells, each cell has to do 1/m job and the death rate is positively correlated with quantity of job. Division rate is independent from job.
p <- 0.05 # p depends on number of living cells
q <- 0.03
weeks <- 1:300
# No. living cell from 1st through 100th week
n.sv.2 <- rep(0,length(weeks))
n.sv.2[1] <- 100
for(i in 2:length(weeks)){
# death rate
p. <- p + 1/n.sv.2[i-1]
# p. <= 1
p. <- min(p.,1)
# death
tmp1 <- sample(0:1,n.sv.2[i-1],replace=TRUE,prob=c(1-p.,p.))
n.death <- sum(tmp1)
# No. living cells for now
tmp2 <- n.sv.2[i-1] - n.death
# Cell division
tmp3 <- sample(0:1,tmp2,replace=TRUE,prob=c(1-q,q))
# No. new born cells
tmp4 <- sum(tmp3)
# No. living cells for now
n.sv.2[i] <- tmp2 + tmp4
}
plot(n.sv.2)
加速なしの場合と並べてプロットしてみる。 Plot two simulation results together.
matplot(cbind(n.sv.1,n.sv.2),type="b")
「今」だけとは限らないが、それでもマルコフ連鎖 MCMC may depend on “NOW” and previous
分裂して生まれたばかりの「赤ん坊細胞」は、すぐには分裂できないとするが、次のタイミングでは分裂可能とする。
Baby cells just born do not divide but they become dividable one period later.
この場合は、生細胞を「赤ん坊」か「そうでないか」の２通りで記録する必要がある。
You have to distinct baby cells and non-baby cells.
p <- 0.05
q <- 0.03
weeks <- 1:300
# No. living cell from 1st through 100th week
# baby and non-baby
n.sv.3.nb <- n.sv.3.bb <- rep(0,length(weeks))
n.sv.3.nb[1] <- 100
for(i in 2:length(weeks)){
# death of non-baby
tmp1.nb <- sample(0:1,n.sv.3.nb[i-1],replace=TRUE,prob=c(1-p,p))
n.death.nb <- sum(tmp1.nb)
# death of baby
tmp1.bb <- sample(0:1,n.sv.3.bb[i-1],replace=TRUE,prob=c(1-p,p))
n.death.bb <- sum(tmp1.bb)
# No. living non-baby cells for now
tmp2 <- n.sv.3.nb[i-1] - n.death.nb
# Cell division
tmp3 <- sample(0:1,tmp2,replace=TRUE,prob=c(1-q,q))
# No. new born cells
tmp4 <- sum(tmp3)
# No. living non-baby cells. Living baby cells are now non-baby
n.sv.3.nb[i] <- tmp2 + n.sv.3.bb[i-1]-n.death.bb
# No. newborns
n.sv.3.bb[i] <- tmp4
}
plot(n.sv.3.nb+n.sv.3.bb)
matplot(cbind(n.sv.3.nb,n.sv.3.bb,n.sv.3.nb+n.sv.3.bb),type="l")
行列を使う Use matrices
ある疾患には、４つの状態があるという。急性軽症・急性重症・慢性・癌化の４つである。この４状態に、健康・死亡の２状態を加えて、全部で６状態を考える。
Assume a disease has 4 statuses, acureLight,acuteSevere,chronic,cancer. Along these 4, two more statuses, healthy and dead, should be considered.
状態に番号を振る。健康：１、急性軽症：２、急性重症：３、慢性：４、癌化：５、死亡：６とする。 Call 4+6 statuses with 1,2,…,6 as 1:healthy, 2:acuteLight, 3:acuteSevere, 4:chronic, 5:cancer, 6:dead.
健康からは、急性軽症、急性重症、死亡(この疾患とは無関係の死亡)に変わりうるとして、その確率をp12,p13,p16とする。
Healthy individual can become acuteLight/acuteSevere/Dead(dead from some reason away from this disease). Let p12,p13,p16 denote their probability.
変わりえないのは、慢性・癌化の２状態であり、その確率p14=p15=0である。 今、健康なまま、という確率はp11とする。
Healthy individual can not get “chronic/cancer” directry, meaning p14=p15=0. (S)he may stay healthy with probability p11.
このとき、p11 + p12 + p13 + p14 + p15 + p16=1であって、p1iはすべて0以上である。これを、状態１からの推移確率ベクトルと呼ぶ。
p = (p11,…,p16) is called the transition vector from status 1 and p11 + … + p16 = 1 and p1i >= 0.
同様にして、1-6のそれぞれに推移確率ベクトルを作ることができる。
Each status has its own transition vector.
6つの状態の状態推移ベクトルを合わせると、6x6行列が作れる。これを状態推移行列という。
With these 6 transition vectors, you can make a 6x6 matrix, called “transition matrix”.
状態推移行列を使って、6状態にある確率の推移を見てみる。
Using this transition matrix, history of probabilities in 6 statuses can be simulated.
推移確率は１年間での状態変化確率とする。
Transition probability is per-year one.
まず行列を作る。Make a matrix.
# From 1:healthy
p1 <- c(0.9,0.08,0.0001,0,0,0.0199)
# From 2:acuteLight
p2 <- c(0.95,0,0,0.04,0,0.01) # 治るか慢性化するか、死亡するか
# From 3:acuteSevere
p3 <- c(0.8,0,0,0.01,0,0.19) # 死亡率0.19
# From 4:chronic
p4 <- c(0.01,0,0,0.95,0.03,0.01) # 癌化する
# From 5:cancer
p5 <- c(0,0,0,0,0.9,0.1)
# From 6: dead
p6 <- c(0,0,0,0,0,1)
# 行列化 matrix
P <- cbind(p1,p2,p3,p4,p5,p6)
P## p1 p2 p3 p4 p5 p6
## [1,] 0.9000 0.95 0.80 0.01 0.0 0
## [2,] 0.0800 0.00 0.00 0.00 0.0 0
## [3,] 0.0001 0.00 0.00 0.00 0.0 0
## [4,] 0.0000 0.04 0.01 0.95 0.0 0
## [5,] 0.0000 0.00 0.00 0.03 0.9 0
## [6,] 0.0199 0.01 0.19 0.01 0.1 1# 全列が足して１になるかを確認 check sum of columns being 1
apply(P,2,sum)## p1 p2 p3 p4 p5 p6
## 1 1 1 1 1 1初期状態は、1:Healthyに居たとする。状態確率ベクトルが(1,0,0,0,0,0)である、という。
Initially you are healthy(1), then your status vector is (1,0,0,0,0,0).
v1 <- c(1,0,0,0,0,0)
# １年後 One-year later
v2 <- P %*% v1
v2## [,1]
## [1,] 0.9000
## [2,] 0.0800
## [3,] 0.0001
## [4,] 0.0000
## [5,] 0.0000
## [6,] 0.0199# ２年後、two-year later
v3 <- P %*% v2
v3## [,1]
## [1,] 0.886080
## [2,] 0.072000
## [3,] 0.000090
## [4,] 0.003201
## [5,] 0.000000
## [6,] 0.038629# x年後、x-year later
x <- 1:50
v.all <- matrix(0,length(x),length(v1))
v.all[1,] <- v1
for(i in 2:length(x)){
v.all[i,] <- P %*% v.all[i-1,]
}
matplot(v.all,type="l")
治療介入してみよう Evaluate therapeutic interventions
慢性状態を治癒に持ち込む確率が、0.01だが、これを0.1に上げる治療A。 Assume a treatment A that increases p41 from 0.01 to 0.1.
PA <- P
PA[1,4] <- 0.1
PA[4,4] <- PA[4,4]-(0.1-0.01)
PA## p1 p2 p3 p4 p5 p6
## [1,] 0.9000 0.95 0.80 0.10 0.0 0
## [2,] 0.0800 0.00 0.00 0.00 0.0 0
## [3,] 0.0001 0.00 0.00 0.00 0.0 0
## [4,] 0.0000 0.04 0.01 0.86 0.0 0
## [5,] 0.0000 0.00 0.00 0.03 0.9 0
## [6,] 0.0199 0.01 0.19 0.01 0.1 1v.all.A <- matrix(0,length(x),length(v1))
v.all.A[1,] <- v1
for(i in 2:length(x)){
v.all.A[i,] <- PA %*% v.all.A[i-1,]
}
# 治療法Aの効果を死亡状態の割合で比較する
matplot(x,cbind(v.all[,6],v.all.A[,6]),type="l")
癌を治す治療Bとして、治癒成功率を0.1にしてみよう。 Assume a treatment B that cures cancer, p51 = 0.1.
PB <- P
PB[1,5] <- 0.1
PB[5,5] <- PB[5,5]-0.1
PB## p1 p2 p3 p4 p5 p6
## [1,] 0.9000 0.95 0.80 0.01 0.1 0
## [2,] 0.0800 0.00 0.00 0.00 0.0 0
## [3,] 0.0001 0.00 0.00 0.00 0.0 0
## [4,] 0.0000 0.04 0.01 0.95 0.0 0
## [5,] 0.0000 0.00 0.00 0.03 0.8 0
## [6,] 0.0199 0.01 0.19 0.01 0.1 1v.all.B <- matrix(0,length(x),length(v1))
v.all.B[1,] <- v1
for(i in 2:length(x)){
v.all.B[i,] <- PB %*% v.all.B[i-1,]
}
# 治療法Aの効果を死亡状態の割合で比較する
matplot(x,cbind(v.all[,6],v.all.A[,6],v.all.B[,6]),type="l")