cocolog:87218726

右三角分布、左三角分布、右台形分布、左台形分布の平均。簡単に計算できるためか、ググッても載っているところがなかった。 (JRF 1594)

JRF 2017年4月14日 (金) Tweet

右三角分布は、a <= x <= b で密度関数が f(x) = 2 * (x - a) / (b - a)^2 で、それ以外の区間では 0 のもの。平均は、(2 * b + a) / 3 になる。

JRF2017/4/146547

左三角分布は、a <= x <= b で密度関数が f(x) = 2 * (b - x) / (b - a)^2 で、それ以外の区間では 0 のもの。平均は、(2 * a + b) / 3 になる。

JRF2017/4/149482

右台形分布は、a <= x <= b で密度関数が f(x) = (c * (b - a) + 2 * (1 - c) * (x - a)) / (b - a)^2 で、それ以外の区間では 0 のもの。平均は、(2 * a + 4 * b - c * (b - a)) / 6 になる。

JRF2017/4/140236

左台形分布は、a <= x <= b で密度関数が f(x) = ((2 - c) * (b - a) - 2 * (1 - c) * (x - a)) / (b - a)^2 で、それ以外の区間では 0 のもの。平均は、(4 * a + 2 * b + c * (b - a)) / 6 になる。

JRF2017/4/145800

分散はまだ計算してない。

JRF2017/4/148664

↓を参考にした。

『計算機シミュレーションのための確率分布乱数生成法』(四辻 哲章 著, プレアデス出版, 2010年 2013年第2版)
https://www.amazon.co.jp/dp/4903814351
http://7net.omni7.jp/detail/1102943434

JRF2017/4/143890

ちなみにこの本の p.74 に上記、左台形分布の式があるが、一ヶ所 + と - の符号が間違っていた。正誤表は↓にあるようだが、そこにはまだ載っていないようだ。

《「計算機シミュレーションのための確率分布乱数生成法」正誤》
http://www.pleiades-publishing.co.jp/download/data/keis.pdf

JRF2017/4/146195

なお、この本によると、乱数発生はそれぞれ次のようになる。

JRF2017/4/144394

……。

右三角分布。

前処理. s = b - a を計算し、値を保持しておく。

Step 1. 区間 [0,1] の一様乱数 U1 と U2 を独立に発生させる。

Step 2. U1 と U2 の値を比較し、Y = max(U1, U2) とする。

Step 3. X = a + s * Y を所望の値とする。

JRF2017/4/145550

……。

左三角分布。

前処理. s = b - a を計算し、値を保持しておく。

Step 1. 区間 [0,1] の一様乱数 U1 と U2 を独立に発生させる。

Step 2. U1 と U2 の値を比較し、Y = min(U1, U2) とする。

Step 3. X = a + s * Y を所望の値とする。

JRF2017/4/147157

……。

右台形分布。

前処理. p = 1 / (1 - c)、q = p * c / 2、r = q^2、s = b - a を計算し、値を保持しておく。

Step 1. 区間 [0,1] の一様乱数 U を発生させる。

Step 2. Y = - q + sqrt(r + p * U) を計算する。

Step 3. X = a + s Y を所望の乱数として計算する。

JRF2017/4/148048

……。

左台形分布。

前処理. p = 1 / (1 - c)、q = p * (2 - c) / 2、r = q^2、s = b - a を計算し、値を保持しておく。

Step 1. 区間 [0,1] の一様乱数 U を発生させる。

Step 2. Y = q - sqrt(r - p * U) を計算する。

Step 3. X = a + s Y を所望の乱数として計算する。

JRF2017/4/143150

……。

なんでこれが必要と思ったかというと、今、シミュレーションを作っていて、ある「率」を平均を指定してランダムに作成したい…つまり、定義域が 0 から 1 の間で少し前までは一様乱数で発生していたのを、平均を動かして使いたいと考えたからである。

JRF2017/4/149991

それで直線的なものとして左三角分布を考えようとしたら、平均は、上で a = 0、b = 1 として 1/3 に決まってしまって好きな平均を出すという目的には合わない。だったら、左台形分布にすれば (2+c)/6、右台形分布にすれば (4-c)/6 が平均になるから、r を平均とすると 0 <= r < 2/3 のとき右台形分布で c = 4-6*r、1/3 < r <= 1 のとき左台形分布で c = 6 * r - 2 にすればいい。範囲が重なっているので r = 0.5 で分ければいいかもしれない。…と考えた。

JRF2017/4/145563

ただ、実際、この乱数を使うかどうかはまだ決めていない。

JRF2017/4/141571

……。

「平均」というより「期待値」と書くべきだったかな…。

JRF2017/4/145294

……。

あっ、ダメだ。

上の論に基づいて台形分布を実際に出力させてみたら、間違いがわかった。まず、平均 r = 0.5 のとき、左台形分布でも右台形分布でも c = 1 となり、上のアルゴリズムのところで p = 1/(1-c) が発散してしまう(ゼロ除算でエラーになる)。r = 0.5 のときは一様分布を使えということだろう。

JRF2017/4/149977

そして、他もいろいろおかしい。例えば、r = 0.2 で右台形分布のアルゴリズムを使うと平均は 1.32... とかおかしな値になるし、試しに左台形分布のアルゴリズムを使ってみると 0.22... とそれっぽい値になるが、gnuplot で分布を描くと x >= 0.6 が 0 の分布になっている。

JRF2017/4/142536

同様に r = 0.8 で左台形分布のアルゴリズムを使うと平均は -0.33... とかおかしな値になるし、試しに右台形分布のアルゴリズムを使ってみると 0.77... とそれっぽい値になるが、gnuplot で分布を描くと x <= 0.4 が 0 の分布になっている。

JRF2017/4/147297

どうも 0 <= c < 1 の範囲しかやってはいけないようだ。よって右台形分布は 0.5 < r <= 2/3、左台形分布は 1/3 <= r < 0.5(、r = 0.5 のとき一様分布)としなければダメなようだ。

JRF2017/4/143930

…そうか、まぁ、今作っているシミュレーションで平均を 1/3 <= r <= 2/3 に収めるというのはまだ修正がきくから、別にいいけど、不便だなぁ…。

JRF2017/4/146228