三角分布に従う乱数

三角分布に従う乱数を発生させるプログラム

最頻値ｇで、0.8ｇ～2.5ｇの範囲の三角分布（当然面積は１です）に従う乱数を発生させるプログラムを
書きたいと思っています。

この三角分布の確率密度関数をＰ（ｘ）とすると、三角分布であるので最頻値ｇの左側である、傾きが正の直線ｈ（ｘ）と
最頻値ｇの右側である、傾きが負の直線ｆ（ｘ）で表せますよね

分布に従う乱数を発生させるためには、これら直線の関数を積分したものの逆関数ｘ＝Ｐ^-1（u） （Ｐのインバースです）
（uは区間[0, 1]の一様乱数）とすればいいというところまでわかったんですが、

とりあえずｈ（ｘ）とｆ（ｘ）をそれぞれ積分して逆関数Ｈ^-1（u）、Ｆ^-1（u）を求めたところまではいいんですが
ｘ＝Ｈ^-1（u）＋Ｆ^-1（u）としてプログラムを実行すると、最頻値ｇの２倍あたりの値（例えば２０に対して３９など）
しか出ず、最頻値ｇより小さい値が出ません。
Ｈ^-1（u）＋Ｆ^-1（u）としているのがダメだと思うのですが、逆関数が２つある場合、ここからどうすればいいですか？

また、初歩的な質問なのですが、区間[0, 1]の一様乱数というのはどう記述すればよいですか？


ぜひ多くの方の回答お待ちしています。
よろしくお願いします。
（最頻値ｇは入力で与えるものとします）

No.2 ベストアンサー 回答者： quaestio 回答日時： 2010/01/31 20:31

最頻値がgで[a, b]の範囲で三角分布に従う乱数を発生させたいとします。

まず、確率密度関数を求めると、
(-∞, 0] : 0
(a, g] : 2(x-a)/{(g-a)(b-a)}
(g, b] : 2(x-b)/{(g-b)(b-a)}
(b, ∞) : 0
分布関数は確率密度関数を積分して、
(-∞, 0] : 0
(a, g] : (x-a)^2/{(g-a)(b-a)}
(g, b] : (g-a)/(b-a)+(x^2-2bx-g^2+2bg)/{(g-b)(b-a)}
(b, ∞) : 1
分布関数の逆関数を求めると、
(-∞, 0] : a
(0, (g-a)/(b-a)] : a+sqrt((g-a)(b-a)x)
((g-a)/(b-a), 1] : b-sqrt((b-g)^2-(b-g)(b-a)(x-(g-a)/(b-a)))
(1, ∞) : b
となりますので、xは0以上1以下の値が入っているとすれば、
if (x <= (g-a)/(b-a)) return a+sqrt((g-a)(b-a)x);
else return b-sqrt((b-g)^2-(b-g)(b-a)(x-(g-a)/(b-a)));
とでも書けばいいでしょう。
（実際に書く場合は、もう少し計算量が少なくなるように書くとは思いますが）

計算があっているかどうかは確認しておいてください。

この回答へのお礼

分かりやすい回答をどうもありがとうございました。
おかげでちゃんと実行することができました。
私のはそもそも積分が間違っていたようです。
詳しい解説、本当にありがとうございました。

お礼日時：2010/02/01 11:00

No.4 回答者： Ishiwara 回答日時： 2010/02/04 10:56

ふつうの長方形分布の乱数を発生させ、条件に合わなければ捨ててしまう、というやり方が一番簡単です。

No.3 回答者： quaestio 回答日時： 2010/01/31 20:40

ANo.2

思いっきり間違いました。
こちらが正しいです。
if (x <= (g-a)/(b-a)) return a+sqrt((g-a)*(b-a)*x);
else return b-sqrt((b-g)*(b-g)-(b-g)*(b-a)*(x-(g-a)/(b-a)));

No.1 回答者： quaestio 回答日時： 2010/01/31 16:40

言語は何を使われているのかわかりませんが、if文等で条件分岐をするのは駄目なのでしょうか？

> また、初歩的な質問なのですが、区間[0, 1]の一様乱数というのはどう記述すればよいですか？

rand関数というような関数がありませんか？
名前は違うかもしれませんが、必ず乱数を発生させる関数があるはずです。

ひょっとして、0からnまでの整数の乱数は得られるけど、[0, 1]の一様乱数をどうやって得るのかわからないということでしょうか？
（C言語なら(double)rand() / RAND_MAXとでもすればいいですが）

この回答へのお礼

回答ありがとうございます。
言語はＣ＋＋です。乱数わかりました。

if文等で条件分岐をするというのは
どういう条件のときにどういう動作をすればよいのかが難しいです・・

お礼日時：2010/01/31 18:40