Applied Statistics and Data Mining

午前の結果に基づいて分布を選択し、パラメータの推定値を与える。

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

# 寿命分布のパラメータ
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.bar <- function(x) {
  return(1-pdensfun(x,rdensfun.params[1],rdensfun.params[2]))
}

# コストレートの式 (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))
}

# nlminbは最小化をしてくれるRの関数
nlminb(1,g)

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

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

# 準備のための乱数データ発生
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))