廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

再生過程における余命の分布(均衡分布)

状況

例えば電球が切れたら交換し切れたら交換し……というようなプロセスを考えます。

電球それぞれの寿命 X_i>0 は独立に同分布 F(x) に従うとします。X_i は非負の連続型確率変数です。

電球の寿命切れ、というイベントが繰り返し生起するため、このような過程を再生過程(renewal process)と呼びます。

このとき、ある時刻 t から観測される電球の余命 W はどのような分布に従うでしょうか。

f:id:abrahamcow:20150418000036p:plain

電球の設置から電球の寿命切れまではすべて同分布 F(x) に従いますが、上図の通り、ある時刻 t から観測される電球の余命は gap time がちょん切られる形になるため、一般には分布 F(x) には従いません。

結論をいうと電球の余命 W は、(1-F(w))/\mu に従います。ここで μ は分布 F(.) の平均です。この分布を近郊分布(equilibrium distribution)と呼びます。

このことの証明は、Cox & Isham (1980) などに与えられています。

この結果を R によるシミュレーションで確かめてみます。

シミュレーション

アルゴリズム

今回行うシミュレーションのアルゴリズムは大雑把にいうと以下の通りです。

  • 乱数 X_i, (i=1,2,\ldots,10000) の累積和をとる
  • 時刻 t は一様乱数によって定める
  • 累積和が t を超える最小の ik とする
  • \sum^{k}_{i=1}X_i - tW の観測値 w として出力
  • 上記を10000回繰り返す

指数分布

ここがパラドキシカルでおもしろいところですが、F(x) が平均 λ の指数分布の場合、(1-F(w))/\mu は平均 λ の指数分布の分布関数と一致します。

設置後一定時間経過しているはずの電球の余命の分布が、設置してからの電球の寿命の分布と等しくなってしまいます。

あらためて考えると、指数分布は無記憶性があるので、このことはさほど驚くべき事実ではないのかもしれません。

逆に (1-F(x))/\mu = F(x) と置き、両辺を微分して方程式を解けば、F (x) =1- \exp(−\lambda x) が得られ、両者が一致するのは X が指数分布に従う場合に限られることがわかります。

#R のコード
n <- 10000
Z_e <- numeric(n)
for(i in 1:n){
  len <- 10000
  x <- rexp(len,1/3)
  t2 <- cumsum(x)
  t1 <- c(0,t2[-len])
  o <- runif(1,0,t2[len])
  sta <-which.max(t2 > o)
  Z_e[i] <-t2[sta] - o #ちょん切る
}
hist(Z_e,freq=FALSE, breaks="scott",xlab="w",main="exponensial")
curve(dexp(x,1/3),add=TRUE)

ヒストグラムと密度関数を重ねてみます。

f:id:abrahamcow:20150417234306p:plain

ガンマ分布

f:id:abrahamcow:20150729175811p:plain

#R のコード
n <- 10000
Z_g <- numeric(n)
for(i in 1:n){
  len <- 10000
  x <- rgamma(len,3,1/2)
  t2 <- cumsum(x)
  t1 <- c(0,t2[-len])
  o <- runif(1,0,t2[len])
  sta <-which.max(t2 > o)
  Z_g[i] <-t2[sta] - o #ちょん切る
}
hist(Z_g,freq=FALSE, breaks="scott",xlab="w",main="gamma")
re_gamma <- function(x){
  pgamma(x,3,1/2, lower.tail=FALSE)/(2*3)
}
curve(re_gamma,add=TRUE,col="red2",lwd=2)

ワイブル分布

#R のコード
Z_W <- numeric(n)
for(i in 1:n){
  len <- 10000
  x <- rweibull(len,2,2)
  t2 <- cumsum(x)
  t1 <- c(0,t2[-len])
  o <- runif(1,0,t2[len])
  sta <-which.max(t2 > o)
  Z_W[i] <-t2[sta] - o
}
hist(Z_W,freq=FALSE,breaks="scott",xlab="w",main="Weibull")
re_wei <- function(x){
  (1-pweibull(x,2,2))/
    (2*gamma(1+1/2))
}
curve(re_wei(x),add=TRUE,col="red2",lwd=2)

f:id:abrahamcow:20150417234525p:plain

対数正規分布

#R のコード
Z_ln <- numeric(n)
for(i in 1:n){
  len <- 10000
  x <- rlnorm(len)
  t2 <- cumsum(x)
  t1 <- c(0,t2[-len])
  o <- runif(1,0,t2[len])
  sta <-which.max(t2 > o)
  Z_ln[i] <-t2[sta] - o
}
hist(Z_ln,freq=FALSE,breaks="scott",xlab="w",main="log-normal")
re_lnorm<- function(x){
  plnorm(x,lower.tail=FALSE)/
    exp(0+1/2)
}
curve(re_lnorm(x),add=TRUE,col="red",lwd=2)

f:id:abrahamcow:20150417234849p:plain

一様分布

#R のコード
Z_u <- numeric(n)
for(i in 1:n){
  len <- 10000
  x <- runif(len,0,1)
  t2 <- cumsum(x)
  t1 <- c(0,t2[-len])
  o <- runif(1,0,t2[len])
  sta <-which.max(t2 > o)
  Z_u[i] <-t2[sta] - o
}
hist(Z_u,freq=FALSE,xlab="w",main="uniform")
re_u <-function(x){
  2-2*x
}
curve(re_u(x),add=TRUE,col="red2",lwd=2)

f:id:abrahamcow:20150417235304p:plain

参考文献

Point Processes (Chapman & Hall/CRC Monographs on Statistics & Applied Probability)

Point Processes (Chapman & Hall/CRC Monographs on Statistics & Applied Probability)

確率過程の基礎

確率過程の基礎

確率モデル入門

確率モデル入門