Applied Statistics and Data Mining

午前の結果に基づいて分布を選択する。 対応する確率分布のpdensfun, ddensfun, rdensfunをRに貼り付ける。

# 寿命分布の累積分布関数
pdensfun = pweibull  # ワイブル分布
pdensfun = plnorm  # 対数正規分布
pdensfun = pgamma  # ガンマ分布

# 寿命分布の確率密度関数
ddensfun = dweibull  # ワイブル分布
ddensfun = dlnorm  # 対数正規分布
ddensfun = dgamma  # ガンマ分布

# 寿命分布の乱数
rdensfun = rweibull  # ワイブル分布
rdensfun = rlnorm  # 対数正規分布
rdensfun = rgamma  # ガンマ分布

パラメータの推定値を与える。

# 寿命分布のパラメータ
rdensfun.params = c(2.2357389,9.7422387)

#今回の想定パラメータ
Cc = 7800+40000
Cp = 7800+5000

時間取替のコストレートの式には、積分が含まれている。Rでは、関数integrateによる積分の機能が提供されている。

三角関数sinを0から2πまで積分してみる。

integrate(sin,0,2*pi)

Rの中で定義済みの関数は、このようにそのまま定積分を計算できる。中で定義されていない関数、あるいは積分変数以外の変数を持つ関数は、積分変数のみを引数として、被積分関数の値を返すRの関数を定義して、integrateに渡す。

f = function(x) {
  f.ret = dlnorm(x, meanlog=2, sdlog=1)
}

このように定義した関数であれば、次のように範囲を指定して、グラフを描くことができる。

plot(f, xlim=c(0,20))

この関数を(0,♾)の範囲で定積分を行うには、次のようにRに指示すれば良い。

integrate(f,0,Inf)

コストレートの式の分子の被積分関数が確率密度関数なので、パラメータまで与えた関数fを定義しておく。

# 確率密度関数
f <- function(x) {
  return(ddensfun(x,rdensfun.params[1],rdensfun.params[2]))
}

コストレートの式の分母の被積分関数が生存関数(=1-累積分布関数のこと)なので、これもパラメータまで与えた関数F.barを定義しておく。

# 生存関数
F.bar <- function(x) {
  return(1-pdensfun(x,rdensfun.params[1],rdensfun.params[2]))
}

コストレートは、分母も分子も積分を含んでいる。上に紹介した関数integrateと、被積分関数の与え方に従えば、次のようにコストレートをRの関数で表現できる。

# コストレートの式 (integrateは1変数関数の数値積分をしてくれるRの関数)
g <- function(x) {
  return((Cc*integrate(f,0,x)$value+Cp*(1-integrate(f,0,x)$value))/(integrate(F.bar,0,x)$value))
}

例えば

g(1)

と入力すると、T=1の場合のコストレートが計算される。

コストレートの最小化に、関数nlminbを使うことができる。nlminbは最小化をしてくれるRの関数で、初期値T=1から始めて、最小値を探すには、次の1行を実行する。

nlminb(1,g)

最小解のみでなく、どのような関数を最小化したかを確認するには、グラフも描いてみると良いだろう。

plot(g,xlim=c(0,20))

他に、次のようにしても、同じグラフを描くことができる。

# コストレートのグラフを描く
g.list <- NULL
for( i in c(1:200) ) {
  g.list <- append(g.list,g(i/10))
}
plot(c(1:500)/10,g.list,type="b")

0.1ずつ関数の値を評価して、折れ線グラフを描いている。

# 再生関数をモンテカルロ近似で求める
# 再生関数が数式で表せる分布の方が少ない

# 準備のための乱数データ発生
M.t.prep <- function() {
  rdist <- function(n) {
    return(rdensfun(n,rdensfun.params[1],rdensfun.params[2]))
  }
  n <- 100000
  X <- rep(0,n)
  K <- 50
  F.k <- array(0,dim=c(K,n))
  for( k in c(1:K) ) {
    X <- X + rdist(n)
    F.k[k,] <- X
  }
  return(F.k)  
}

# 再生関数をモンテカルロ計算で求める
M.t <- function(block.cycles) {
  n <- 100000
  X <- rep(0,n)
  K <- 50
  M.ret <- rep(0,length(block.cycles))
  for( i in c(1:length(block.cycles)) ) {
    for( k in c(1:K) ) {
      M.ret[i] <- M.ret[i] + sum(M.t.data[k,]<=block.cycles[i])/n
    }
  }
  return(M.ret)
}

# 使う前に一度だけ準備を行う（分布を変更する都度、実行しなおす必要あり）
M.t.data = M.t.prep()

# 再生関数のグラフ
plot(1:100,M.t(1:100))

# ブロック取り替えのコストレートのグラフ
plot(1:100,(M.t(1:100)*Cc+Cp)/(1:100),
     xlab="T",
     ylab="Cost Rate")

# ブロック取り替えの子ストレート
Cost.rate.BR <- function(cycle.length) {
  return((M.t(cycle.length)*Cc+Cp)/cycle.length)
}

# ブロック取り替えの最適化をRの関数optimizeで実行する
optimize(Cost.rate.BR,c(0,10))