差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
r:maintenance:hardtime [2018/12/10 08:46] – watalur:maintenance:hardtime [2018/12/10 09:18] (現在) – [時間計画保全] watalu
行 5: 行 5:
  
 <code> <code>
-select.distribution <- function(dist="weibull") {+select.distribution = function(dist="weibull") {
   if(dist %in% c("weibull","gamma","lnorm")) {   if(dist %in% c("weibull","gamma","lnorm")) {
     distrname <- dist     distrname <- dist
行 29: 行 29:
  
 <code> <code>
-model <- select.distribution("lognormal") +model = select.distribution("lognormal") 
-model <- select.distribution("weibull") +model = select.distribution("weibull") 
-model <- select.distribution("gamma")+model = select.distribution("gamma")
 </code> </code>
  
行 87: 行 87:
 <code> <code>
 # 確率密度関数 # 確率密度関数
-f <- function(x) {+f = function(x) {
   return(model$d(x,model.parameters[1],model.parameters[2]))   return(model$d(x,model.parameters[1],model.parameters[2]))
 } }
行 95: 行 95:
 <code> <code>
 # 生存関数 # 生存関数
-F.bar <- function(x) {+F.bar = function(x) {
   return(1-model$p(x,model.parameters[1],model.parameters[2]))   return(1-model$p(x,model.parameters[1],model.parameters[2]))
 } }
行 106: 行 106:
 <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))
 } }
行 134: 行 134:
 <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")
行 153: 行 153:
 <code> <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
行 173: 行 173:
 <code> <code>
 # 再生関数をモンテカルロ計算で求める # 再生関数をモンテカルロ計算で求める
-M.t <- function(block.cycles) {+M.t = function(block.cycles) {
   n <- 100000   n <- 100000
   X <- rep(0,n)   X <- rep(0,n)
行 192: 行 192:
 # 使う前に一度だけ準備を行う（分布を変更する都度、実行しなおす必要あり） # 使う前に一度だけ準備を行う（分布を変更する都度、実行しなおす必要あり）
 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>