差分

このページの2つのバージョン間の差分を表示します。

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
r:maintenance:hardtime [2018/12/10 07:40] – [生存関数] watalur:maintenance:hardtime [2018/12/10 09:18] (現在) – [時間計画保全] watalu
行 2: 行 2:
 === 寿命分布の設定 === === 寿命分布の設定 ===
  
-午前の結果に基づいて分布を選択し、パラメータの推定値を与える。+午前の結果に基づいて分布を選択する関数をselect.distributionとして与えておく。 
 <code> <code>
-# 寿命分布 +select.distribution = function(dist="weibull") { 
-pdensfun = pweibull  # ワイブル分布 +  if(dist %in% c("weibull","gamma","lnorm")) { 
-pdensfun = plnorm  # 対数正規分布 +    distrname <- dist 
-pdensfun = pgamma  # ガンマ分布+    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; 
 +  } 
 +} 
 +</code>
  
 +対数正規分布、ワイブル分布、ガンマ分布のいずれかを選択するために、次の３行から１行を選んで実行する。
 +
 +<code>
 +model = select.distribution("lognormal")
 +model = select.distribution("weibull")
 +model = select.distribution("gamma")
 +</code>
 +
 +同様に、午前に推定したパラメータの推定値を与える。
 +
 +<code>
 # 寿命分布のパラメータ # 寿命分布のパラメータ
-rdensfun.params = c(2.2357389,9.7422387)+model.parameters = c(2.2357389,9.7422387)
 </code> </code>
 +
 +これら２つのオブジェクトmodel、model.parametersは後で用いるので、名称を変更したら、下のコードも変更が必要となる。
  
 === コストの設定 === === コストの設定 ===
行 21: 行 51:
 </code> </code>
  
-==== 生存関数 ===+これらもオブジェクト名は変更しない方がいい。
  
 +==== 数値積分 ====
 +
 +時間取替のコストレートの式には、積分が含まれている。Rでは、関数integrateによる積分の機能が提供されている。
 +
 +三角関数sinを0から2πまで積分してみる。
 +
 +<code>
 +integrate(sin,0,2*pi)
 +</code>
 +
 +Rの中で定義済みの関数は、このようにそのまま定積分を計算できる。中で定義されていない関数、あるいは積分変数以外の変数を持つ関数は、積分変数のみを引数として、被積分関数の値を返すRの関数を定義して、integrateに渡す。
 +
 +<code>
 +f = function(x) {
 +  f.ret = dlnorm(x, meanlog=2, sdlog=1)
 +  return(f.ret)
 +}
 +</code>
 +
 +このように定義した関数であれば、次のように範囲を指定して、グラフを描くことができる。
 +<code>
 +plot(f, xlim=c(0,20))
 +</code>
 +
 +この関数を(0,♾)の範囲で定積分を行うには、次のようにRに指示すれば良い。
 +<code>
 +integrate(f,0,Inf)
 +</code>
 +
 +==== 確率分布の密度関数と生存関数 ====
 +
 +コストレートの式の分子の被積分関数が確率密度関数なので、パラメータまで与えた関数fを定義しておく。
 +<code>
 +# 確率密度関数
 +f = function(x) {
 +  return(model$d(x,model.parameters[1],model.parameters[2]))
 +}
 +</code>
 +
 +コストレートの式の分母の被積分関数が生存関数(=1-累積分布関数のこと)なので、これもパラメータまで与えた関数F.barを定義しておく。
 <code> <code>
 # 生存関数 # 生存関数
-F.bar <- function(x) { +F.bar = function(x) { 
-  return(1-pdensfun(x,rdensfun.params[1],rdensfun.params[2]))+  return(1-model$p(x,model.parameters[1],model.parameters[2]))
 } }
 </code> </code>
  
-=== 時間取り替え ===+==== 時間取り替え ==== 
 + 
 +コストレートは、分母も分子も積分を含んでいる。上に紹介した関数integrateと、被積分関数の与え方に従えば、次のようにコストレートをRの関数で表現できる。
  
 <code> <code>
 # コストレートの式 (integrateは1変数関数の数値積分をしてくれるRの関数) # コストレートの式 (integrateは1変数関数の数値積分をしてくれるRの関数)
-g <- function(x) {+g = function(x) {
   return((Cc*integrate(f,0,x)$value+Cp*(1-integrate(f,0,x)$value))/(integrate(F.bar,0,x)$value))   return((Cc*integrate(f,0,x)$value+Cp*(1-integrate(f,0,x)$value))/(integrate(F.bar,0,x)$value))
 } }
 </code> </code>
  
 +例えば
 +<code>
 +g(1)
 +</code>
 +と入力すると、T=1の場合のコストレートが計算される。
 +
 +コストレートの最小化に、関数nlminbを使うことができる。nlminbは最小化をしてくれるRの関数で、初期値T=1から始めて、最小値を探すには、次の1行を実行する。
 <code> <code>
-# nlminbは最小化をしてくれるRの関数 
 nlminb(1,g) nlminb(1,g)
 +</code>
 +関数optimizeを使うこともできる。こちらは、探索する変数の範囲を最小値と最大値で指定する。以下ではc(0,10)として、0以上10以下で探索させてみた。
 +<code>
 +optimize(g,c(0,10))
 </code> </code>
  
 +最小解のみでなく、どのような関数を最小化したかを確認するには、グラフも描いてみると良いだろう。
 +<code>
 +plot(g,xlim=c(0,20))
 +</code>
 +これは、関数gがplotで描けない流儀で定義されているため、エラーになる。
 +この関数のグラフを描くには、次のようにする。
 <code> <code>
 # コストレートのグラフを描く # コストレートのグラフを描く
-g.list <- NULL+g.list = NULL
 for( i in c(1:500) ) { for( i in c(1:500) ) {
-  g.list <- append(g.list,g(i/10))+  g.list = append(g.list,g(i/10))
 } }
 plot(c(1:500)/10,g.list,type="b") plot(c(1:500)/10,g.list,type="b")
 </code> </code>
 +0.1から50.0まで、0.1ずつ関数の値を評価して、折れ線グラフを描いている。
  
-=== ブロック取り替え ===+==== ブロック取り替え ====
  
-<code> +ブロック取替のコストレートの式には、再生過程の再生関数が含まれている。これは再生分布が正規分布か指数分布の場合ぐらいしか、解析的に表現することができない。
-# 再生関数をモンテカルロ近似で求める +
-# 再生関数が数式で表せる分布の方が少ない+
  
 +そこで、モンテカルロ法を用いて、近似することにする。モンテカルロ法は、モンテカルロ積分とも呼ばれる。定積分を、確率分布に関する期待値の計算と解釈できるように変形する。そして、その分布に従う乱数を発生させて代入した、被積分関数の値をたくさん計算し、それらの平均値を定積分の近似値とする。
 +
 +まずはモンテカルロ積分に用いる乱数をたくさん生成する関数M.t.prepを定義する。
 +ここで生成するのは、１回目の再生までの時間の乱数、２回目の再生までの時間の乱数、などであり、５０回目まで用意しておくように、コードを準備してある。
 +
 +<code>
 # 準備のための乱数データ発生 # 準備のための乱数データ発生
-M.t.prep <- function() {+M.t.prep = function() {
   rdist <- function(n) {   rdist <- function(n) {
-    return(rdensfun(n,rdensfun.params[1],rdensfun.params[2]))+    return(model$r(n,model.parameters[1],model.parameters[2]))
   }   }
   n <- 100000   n <- 100000
行 74: 行 167:
   return(F.k)     return(F.k)  
 } }
 +</code>
  
 +この関数の生成した乱数をオブジェクトM.t.dataに保管する前提で、次のモンテカルロ積分による関数M.tは再生関数の計算を行う。
 +
 +<code>
 # 再生関数をモンテカルロ計算で求める # 再生関数をモンテカルロ計算で求める
-M.t <- function(block.cycles) {+M.t = function(block.cycles) {
   n <- 100000   n <- 100000
   X <- rep(0,n)   X <- rep(0,n)
行 89: 行 186:
 } }
 </code> </code>
 +
 +上の２つの関数を用いて、再生関数を計算したり、ブロック取替のコストレートを計算したり、それを最適化するためのＲのコードは以下の通りになる。
  
 <code> <code>
 # 使う前に一度だけ準備を行う（分布を変更する都度、実行しなおす必要あり） # 使う前に一度だけ準備を行う（分布を変更する都度、実行しなおす必要あり）
 M.t.data = M.t.prep() 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:100,M.t(1:100))+plot((1:300)/10,M.t((1:300)/10), 
 +     xlab="Time",  
 +     ylab="M(t)")
  
 # ブロック取り替えのコストレートのグラフ # ブロック取り替えのコストレートのグラフ
-plot(1:100,(M.t(1:100)*Cc+Cp)/(1:100),+plot((1:300)/10,(M.t((1:300)/10)*Cc+Cp)/(1:300)/10,
      xlab="T",      xlab="T",
      ylab="Cost Rate")      ylab="Cost Rate")
  
 # ブロック取り替えの子ストレート # ブロック取り替えの子ストレート
-Cost.rate.BR <- function(cycle.length) {+g = function(cycle.length) {
   return((M.t(cycle.length)*Cc+Cp)/cycle.length)   return((M.t(cycle.length)*Cc+Cp)/cycle.length)
 } }
  
 # ブロック取り替えの最適化をRの関数optimizeで実行する # ブロック取り替えの最適化をRの関数optimizeで実行する
-optimize(Cost.rate.BR,c(0,10))+optimize(g,c(0,10))
 </code> </code>