差分
このページの2つのバージョン間の差分を表示します。
両方とも前のリビジョン前のリビジョン次のリビジョン | 前のリビジョン | ||
r:maintenance:hardtime [2018/12/10 07:53] – [数値積分] watalu | r:maintenance:hardtime [2018/12/10 09:18] (現在) – [時間計画保全] watalu | ||
---|---|---|---|
行 2: | 行 2: | ||
=== 寿命分布の設定 === | === 寿命分布の設定 === | ||
- | 午前の結果に基づいて分布を選択し、パラメータの推定値を与える。 | + | 午前の結果に基づいて分布を選択する関数をselect.distributionとして与えておく。 |
< | < | ||
- | # 寿命分布 | + | select.distribution |
- | pdensfun | + | if(dist %in% c(" |
- | pdensfun | + | |
- | pdensfun | + | pdistr <- eval(parse(text=paste(" |
+ | rdistr <- eval(parse(text=paste(" | ||
+ | ddistr <- eval(parse(text=paste(" | ||
+ | qdistr <- eval(parse(text=paste(" | ||
+ | return(list(p=pdistr, | ||
+ | } else if (dist == " | ||
+ | | ||
+ | pdistr <- eval(parse(text=paste(" | ||
+ | rdistr <- eval(parse(text=paste(" | ||
+ | ddistr <- eval(parse(text=paste(" | ||
+ | qdistr <- eval(parse(text=paste(" | ||
+ | return(list(p=pdistr, | ||
+ | } else { | ||
+ | break; | ||
+ | } | ||
+ | } | ||
+ | </ | ||
+ | 対数正規分布、ワイブル分布、ガンマ分布のいずれかを選択するために、次の3行から1行を選んで実行する。 | ||
+ | |||
+ | < | ||
+ | model = select.distribution(" | ||
+ | model = select.distribution(" | ||
+ | model = select.distribution(" | ||
+ | </ | ||
+ | |||
+ | 同様に、午前に推定したパラメータの推定値を与える。 | ||
+ | |||
+ | < | ||
# 寿命分布のパラメータ | # 寿命分布のパラメータ | ||
- | rdensfun.params | + | model.parameters |
</ | </ | ||
+ | |||
+ | これら2つのオブジェクトmodel、model.parametersは後で用いるので、名称を変更したら、下のコードも変更が必要となる。 | ||
=== コストの設定 === | === コストの設定 === | ||
行 20: | 行 50: | ||
Cp = 7800+5000 | Cp = 7800+5000 | ||
</ | </ | ||
+ | |||
+ | これらもオブジェクト名は変更しない方がいい。 | ||
==== 数値積分 ==== | ==== 数値積分 ==== | ||
行 36: | 行 68: | ||
f = function(x) { | f = function(x) { | ||
f.ret = dlnorm(x, meanlog=2, sdlog=1) | f.ret = dlnorm(x, meanlog=2, sdlog=1) | ||
+ | return(f.ret) | ||
} | } | ||
</ | </ | ||
行 49: | 行 82: | ||
</ | </ | ||
- | ==== 生存関数 === | + | ==== 確率分布の密度関数と生存関数 |
+ | コストレートの式の分子の被積分関数が確率密度関数なので、パラメータまで与えた関数fを定義しておく。 | ||
+ | < | ||
+ | # 確率密度関数 | ||
+ | f = function(x) { | ||
+ | return(model$d(x, | ||
+ | } | ||
+ | </ | ||
+ | |||
+ | コストレートの式の分母の被積分関数が生存関数(=1-累積分布関数のこと)なので、これもパラメータまで与えた関数F.barを定義しておく。 | ||
< | < | ||
# 生存関数 | # 生存関数 | ||
- | 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])) |
} | } | ||
</ | </ | ||
+ | |||
==== 時間取り替え ==== | ==== 時間取り替え ==== | ||
+ | |||
+ | コストレートは、分母も分子も積分を含んでいる。上に紹介した関数integrateと、被積分関数の与え方に従えば、次のようにコストレートをRの関数で表現できる。 | ||
< | < | ||
# コストレートの式 (integrateは1変数関数の数値積分をしてくれるRの関数) | # コストレートの式 (integrateは1変数関数の数値積分をしてくれるRの関数) | ||
- | g <- function(x) { | + | g = function(x) { |
return((Cc*integrate(f, | return((Cc*integrate(f, | ||
} | } | ||
</ | </ | ||
+ | 例えば | ||
+ | < | ||
+ | g(1) | ||
+ | </ | ||
+ | と入力すると、T=1の場合のコストレートが計算される。 | ||
+ | |||
+ | コストレートの最小化に、関数nlminbを使うことができる。nlminbは最小化をしてくれるRの関数で、初期値T=1から始めて、最小値を探すには、次の1行を実行する。 | ||
< | < | ||
- | # nlminbは最小化をしてくれるRの関数 | ||
nlminb(1,g) | nlminb(1,g) | ||
+ | </ | ||
+ | 関数optimizeを使うこともできる。こちらは、探索する変数の範囲を最小値と最大値で指定する。以下ではc(0, | ||
+ | < | ||
+ | optimize(g, | ||
</ | </ | ||
+ | 最小解のみでなく、どのような関数を最小化したかを確認するには、グラフも描いてみると良いだろう。 | ||
+ | < | ||
+ | plot(g, | ||
+ | </ | ||
+ | これは、関数gがplotで描けない流儀で定義されているため、エラーになる。 | ||
+ | この関数のグラフを描くには、次のようにする。 | ||
< | < | ||
# コストレートのグラフを描く | # コストレートのグラフを描く | ||
- | g.list | + | g.list |
for( i in c(1:500) ) { | for( i in c(1:500) ) { | ||
- | g.list | + | g.list |
} | } | ||
plot(c(1: | plot(c(1: | ||
</ | </ | ||
+ | 0.1から50.0まで、0.1ずつ関数の値を評価して、折れ線グラフを描いている。 | ||
==== ブロック取り替え ==== | ==== ブロック取り替え ==== | ||
- | < | + | ブロック取替のコストレートの式には、再生過程の再生関数が含まれている。これは再生分布が正規分布か指数分布の場合ぐらいしか、解析的に表現することができない。 |
- | # 再生関数をモンテカルロ近似で求める | + | |
- | # 再生関数が数式で表せる分布の方が少ない | + | |
+ | そこで、モンテカルロ法を用いて、近似することにする。モンテカルロ法は、モンテカルロ積分とも呼ばれる。定積分を、確率分布に関する期待値の計算と解釈できるように変形する。そして、その分布に従う乱数を発生させて代入した、被積分関数の値をたくさん計算し、それらの平均値を定積分の近似値とする。 | ||
+ | |||
+ | まずはモンテカルロ積分に用いる乱数をたくさん生成する関数M.t.prepを定義する。 | ||
+ | ここで生成するのは、1回目の再生までの時間の乱数、2回目の再生までの時間の乱数、などであり、50回目まで用意しておくように、コードを準備してある。 | ||
+ | |||
+ | < | ||
# 準備のための乱数データ発生 | # 準備のための乱数データ発生 | ||
- | M.t.prep | + | M.t.prep |
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 | ||
行 101: | 行 167: | ||
return(F.k) | return(F.k) | ||
} | } | ||
+ | </ | ||
+ | この関数の生成した乱数をオブジェクトM.t.dataに保管する前提で、次のモンテカルロ積分による関数M.tは再生関数の計算を行う。 | ||
+ | |||
+ | < | ||
# 再生関数をモンテカルロ計算で求める | # 再生関数をモンテカルロ計算で求める | ||
- | M.t <- function(block.cycles) { | + | M.t = function(block.cycles) { |
n <- 100000 | n <- 100000 | ||
X <- rep(0,n) | X <- rep(0,n) | ||
行 116: | 行 186: | ||
} | } | ||
</ | </ | ||
+ | |||
+ | 上の2つの関数を用いて、再生関数を計算したり、ブロック取替のコストレートを計算したり、それを最適化するためのRのコードは以下の通りになる。 | ||
< | < | ||
# 使う前に一度だけ準備を行う(分布を変更する都度、実行しなおす必要あり) | # 使う前に一度だけ準備を行う(分布を変更する都度、実行しなおす必要あり) | ||
M.t.data = M.t.prep() | M.t.data = M.t.prep() | ||
+ | |||
+ | # 生成した乱数のヒストグラムを10回目ぐらいまで描いてみる | ||
+ | hist(M.t.data[1, | ||
+ | hist(M.t.data[2, | ||
+ | hist(M.t.data[3, | ||
+ | hist(M.t.data[4, | ||
+ | hist(M.t.data[5, | ||
+ | hist(M.t.data[6, | ||
+ | hist(M.t.data[7, | ||
+ | hist(M.t.data[8, | ||
+ | hist(M.t.data[9, | ||
+ | hist(M.t.data[10, | ||
# 再生関数のグラフ | # 再生関数のグラフ | ||
- | plot(1:100,M.t(1:100)) | + | plot((1:300)/10,M.t((1:300)/10), |
+ | | ||
+ | | ||
# ブロック取り替えのコストレートのグラフ | # ブロック取り替えのコストレートのグラフ | ||
- | plot(1:100,(M.t(1:100)*Cc+Cp)/ | + | plot((1:300)/10,(M.t((1:300)/10)*Cc+Cp)/ |
| | ||
| | ||
# ブロック取り替えの子ストレート | # ブロック取り替えの子ストレート | ||
- | Cost.rate.BR <- function(cycle.length) { | + | g = function(cycle.length) { |
return((M.t(cycle.length)*Cc+Cp)/ | return((M.t(cycle.length)*Cc+Cp)/ | ||
} | } | ||
# ブロック取り替えの最適化をRの関数optimizeで実行する | # ブロック取り替えの最適化をRの関数optimizeで実行する | ||
- | optimize(Cost.rate.BR,c(0,10)) | + | optimize(g,c(0,10)) |
</ | </ |