廿TT

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

レイリー分布を仮定した可用性(2 状態システムの稼働率)の推定

状況

なんらかの機械を設置したとする.

この機械はある程度時間が経過すると故障する.

修理にはまたある程度の時間がかかる.

つまりこの機械は, 稼働中と修理中の 2 つの状態を交互に繰り返す.

f:id:abrahamcow:20150107212609p:plain

いま, この機械をある時刻に観測したとき, 機械が稼働している割合に興味がある.

このような状況は交代再生過程(alternating renewal process)によってモデル化され, 機械が稼働している割合を可用性(availability; アベイラビリティ)と呼ぶ.

可用性の推定

故障間隔を確率変数 X_i \overset{iid}{\sim} F_X, 修復時間を確率変数  Y_i \overset{iid}{\sim} F_Y で表す(i=1,\ldots,n).

可用性の極限は,

\displaystyle A = \frac{E[X_i]}{E[X_i]+E[Y_i]}

A の自明な推定量は,

\displaystyle \hat A = \frac{\sum X_i/n}{\sum X_i/n+\sum Y_i/n}=\frac{\sum X_i}{\sum X_i + \sum Y_i}

である.

一方, 故障間隔と修復時間に特定の分布を仮定して, パラメトリックに可用性を推定する方法も考えられる.

シミュレーション

今回は故障間隔, 修復時間にそれぞれパラメータ 8, 2 のレイリー分布(シェイプパラメータ = 2 のワイブル分布)を仮定する.

(本当はワイブル分布を使いたいが, 最尤推定量が閉じた形で求まらないため, 少サンプルだといい結果がでなかった.)

可用性は

 \displaystyle \frac{8}{8+2}=0.8.

上で「自明な推定量」と書いた \hat A による推定を手法 1 と呼ぶことにする.

また故障間隔と修復時間の分布パラメータを \eta_1,\eta_2 として, その最尤推定量を \hat \eta _1,\hat \eta _2 と表記すると,

 \displaystyle \frac{\hat \eta _1}{ \hat \eta _1 + \hat \eta _2}

より可用性を推定できる. これを手法 2 とする.

レイリー分布のパラメータの最尤推定量は,

 \displaystyle \hat \eta =\sqrt{\frac{1}{n}\sum_{i=1}^{n}X_i^2}.

10000回のシミュレーションの結果, 手法 2 のほうが若干 MSE(平均二乗誤差)が小さい模様.

f:id:abrahamcow:20150619201134p:plain

推定値のばらつきも, 手法 2 のほうが少しだけ小さい.

f:id:abrahamcow:20150619201316p:plain

シミュレーションの詳細は以下の R のコードを参照.

eta1 <- 8
eta2 <- 2
shp <- 2
rho <- eta1/(eta1 + eta2)
est1 <- est2 <-numeric(1000)
mse1 <- mse2  <-numeric(1000)
library(MASS)
simu <- function(n){
  z <- numeric(n)
  z[2*(1:(n/2))-1] <- rweibull(n/2,shape=shp,scale=eta1) 
  z[2*(1:(n/2))] <- rweibull(n/2,shape=shp,scale=eta2)
  s <- 1:n %% 2
  data.frame(z,s)
}
n=10
for(i in 1:1000){
  dat1 <-simu(n)
  est1[i] <-sum(dat1$z *dat1$s) / sum(dat1$z)
  mse1[i] <- (est1[i] - rho)^2
  x <- with(dat1,z[s==1])
  y <- with(dat1,z[s==0])
  eta_hat1 <- 
    (sum(x^shp) / (n/2))^(1/shp)
  eta_hat2 <- 
    (sum(y^shp) / (n/2))^(1/shp)
  est2[i] <-eta_hat1 / (eta_hat1 + eta_hat2)
  mse2[i] <- (est2[i] - rho)^2
}
mlim <-max(mse1,mse2)
plot(mse1,mse2,xlim=c(0,mlim),ylim=c(0,mlim))
abline(0,1,col="red3",lwd=2)
boxplot(est1,est2)

参考文献

確率過程の基礎

確率過程の基礎


可用性 - Wikipedia