Applied Statistics and Data Mining

午前の結果に基づいて分布を選択する関数をselect.distributionとして与えておく。

select.distribution = function(dist="weibull") {
  if(dist %in% c("weibull","gamma","lnorm")) {
    distrname <- dist
    pdistr <- eval(parse(text=paste("p",distrname,sep="")))
    rdistr <- eval(parse(text=paste("r",distrname,sep="")))
    ddistr <- eval(parse(text=paste("d",distrname,sep="")))
    qdistr <- eval(parse(text=paste("q",distrname,sep="")))
    return(list(p=pdistr,r=rdistr,d=ddistr,q=qdistr))
  } else if (dist == "lognormal") {
    distrname <- "lnorm"
    pdistr <- eval(parse(text=paste("p",distrname,sep="")))
    rdistr <- eval(parse(text=paste("r",distrname,sep="")))
    ddistr <- eval(parse(text=paste("d",distrname,sep="")))
    qdistr <- eval(parse(text=paste("q",distrname,sep="")))
    return(list(p=pdistr,d=ddistr,q=qdistr,r=rdistr))
  } else {
    break;
  }
}

対数正規分布、ワイブル分布、ガンマ分布のいずれかを選択するために、次の３行から１行を選んで実行する。

model = select.distribution("lognormal")
model = select.distribution("weibull")
model = select.distribution("gamma")

同様に、午前に推定したパラメータの推定値を与える。

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

これら２つのオブジェクトmodel、model.parametersは後で用いるので、名称を変更したら、下のコードも変更が必要となる。

#今回の想定パラメータ
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)
  return(f.ret)
}

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

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

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

integrate(f,0,Inf)

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

# 確率密度関数
f = function(x) {
  return(model$d(x,model.parameters[1],model.parameters[2]))
}

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

# 生存関数
F.bar = function(x) {
  return(1-model$p(x,model.parameters[1],model.parameters[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)

関数optimizeを使うこともできる。こちらは、探索する変数の範囲を最小値と最大値で指定する。以下ではc(0,10)として、0以上10以下で探索させてみた。

optimize(g,c(0,10))

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

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

これは、関数gがplotで描けない流儀で定義されているため、エラーになる。 この関数のグラフを描くには、次のようにする。

# コストレートのグラフを描く
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")

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

ブロック取替のコストレートの式には、再生過程の再生関数が含まれている。これは再生分布が正規分布か指数分布の場合ぐらいしか、解析的に表現することができない。

そこで、モンテカルロ法を用いて、近似することにする。モンテカルロ法は、モンテカルロ積分とも呼ばれる。定積分を、確率分布に関する期待値の計算と解釈できるように変形する。そして、その分布に従う乱数を発生させて代入した、被積分関数の値をたくさん計算し、それらの平均値を定積分の近似値とする。

まずはモンテカルロ積分に用いる乱数をたくさん生成する関数M.t.prepを定義する。 ここで生成するのは、１回目の再生までの時間の乱数、２回目の再生までの時間の乱数、などであり、５０回目まで用意しておくように、コードを準備してある。

# 準備のための乱数データ発生
M.t.prep = function() {
  rdist <- function(n) {
    return(model$r(n,model.parameters[1],model.parameters[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.dataに保管する前提で、次のモンテカルロ積分による関数M.tは再生関数の計算を行う。

# 再生関数をモンテカルロ計算で求める
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()

# 生成した乱数のヒストグラムを10回目ぐらいまで描いてみる
hist(M.t.data[1,],prob=TRUE,breaks=c(0:ceiling(max(M.t.data[1:10,]))),xlim=c(0,max(M.t.data[1:10,])),xlab="Time",main="Renewal Distributions")
hist(M.t.data[2,],prob=TRUE,breaks=c(0:ceiling(max(M.t.data[1:10,]))),add=TRUE)
hist(M.t.data[3,],prob=TRUE,breaks=c(0:ceiling(max(M.t.data[1:10,]))),add=TRUE)
hist(M.t.data[4,],prob=TRUE,breaks=c(0:ceiling(max(M.t.data[1:10,]))),add=TRUE)
hist(M.t.data[5,],prob=TRUE,breaks=c(0:ceiling(max(M.t.data[1:10,]))),add=TRUE)
hist(M.t.data[6,],prob=TRUE,breaks=c(0:ceiling(max(M.t.data[1:10,]))),add=TRUE)
hist(M.t.data[7,],prob=TRUE,breaks=c(0:ceiling(max(M.t.data[1:10,]))),add=TRUE)
hist(M.t.data[8,],prob=TRUE,breaks=c(0:ceiling(max(M.t.data[1:10,]))),add=TRUE)
hist(M.t.data[9,],prob=TRUE,breaks=c(0:ceiling(max(M.t.data[1:10,]))),add=TRUE)
hist(M.t.data[10,],prob=TRUE,breaks=c(0:ceiling(max(M.t.data[1:10,]))),add=TRUE)

# 再生関数のグラフ
plot((1:300)/10,M.t((1:300)/10),
     xlab="Time", 
     ylab="M(t)")

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

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

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