SCILABの正規分布乱数を使ったモンテカルロシミュレーション

作成:2000-04-06

更新:2000-10-03


★正規分布の足し算

正規分布っていうのは←こんなやつで、

←こんなやつですが、

普通の乱数っていうのは0〜1まで均等に分布しますが、工学のシミュレーションをするときには、だいたい、モノの精度って正規分布してますから、正規分布な乱数が欲しいですね。

技術評論社:奥村晴彦著:C言語によるアルゴリズム辞典 によれば rnd()を0〜1.00の一様乱数として、

t=sqrt(-2.0 * log(1-rnd()));
u=2*PI*rnd();

r1=t*cos(u) ; r2=t*sin(u);

で2個ずつ、平均0、分散1の正規分布する乱数が得られるそうです。

SCILABも(多分)これによる正規分布乱数を発生する関数 rand(,,"normal")っていうのが実装されてますので早速使ってみましょう。


//平均0、標準偏差1の正規分布乱数は randn() で求められる。 この乱数を10000個発生させて、

ns = rand(1,10000 , "normal"); // 1×10000 個の行列=1次元配列
//各統計データを求めると

mean(ns) //平均
ans = 0.00132466

st_deviation(ns) 標準偏差
ans = 1.00089

min(ns) 最小値
ans = -3.83109

max(ns) 最大値
ans = 3.84904

ヒストグラムを書いてみよう。
って簡単にヒストグラムが書けるところが便利っすね

histplot([-4:0.1:4] , ns ,5) ; //[-4:0.1:4]がヒストグラムの幅を示す数列 5はカラー(赤)

//標準正規分布じゃなくて、平均 xA 標準偏差 sA の場合は

xA = 3 ; sA = 2 ; nsA = randn(1,10000) * sA + xA;

となります。* sA で各要素を sA倍するんですが、
各要素に xA を足す場合は + xA であって、 .+ xA というMaTXとの記述とは違う事に注意(この辺、ツールで異なるのはしかたがないんだが、なんとかならんもんか?)

この場合、

mean(ns) //平均
ans = 2.97621

st_deviation(ns) //標準偏差
ans = 2.00618

min(ns) //最小値
ans = -4.0327

max(ns) //最大値
ans = 12.5495


んじゃ、応用問題で、

長さ 3cm の棒と 長さ 5cm の棒を足す計算をしてみる。って、
それは 3+5=8cm っすよね。

じゃ、その3cmっていうのが「平均3cm,標準偏差0.1cm
   その5cmっていうのが「平均5cm,標準偏差0.3cm

だとどうだろう?

平均=3.0 標準偏差=0.1 の正規分布と 平均=5.0 標準偏差=0.3 の正規分布とを足して新しい分布を作り、 その分布の平均と標準偏差をシミュレートしてみる。

xA = 3 ; sA = 0.1 ; nsA = rand(1,10000,"normal") * sA + xA;
xB = 5 ; sB = 0.3 ; nsB = rand(1,10000,"normal") * sB + xB;

{ mean(nsA) ,st_deviation(nsA), min(nsA), max(nsA) }
ans = {2.99965, 0.100441, 2.61988, 3.36458}

{ mean(nsB) ,st_deviation(nsB),min(nsB), max(nsB) }
ans = {4.99495, 0.30041 , 3.82579, 6.27571}

//この nsA と nsB から nsC を作る。

nsC = nsA + nsB;

{ mean(nsC) ,std(nsC), min(nsC), max(nsC) }
ans = {7.9946, 0.317144, 6.84426 , 9.25621}

 

{ mean(nsA)   ,std(nsA),  min(nsA), max(nsA) }
//     平均    標準偏差    最小    最大
 ans = {2.99965, 0.100441, 2.61988, 3.36458}

//理論値 3      0.1

 

histplot([1.5:0.1:4.5],nsA,5)

{ mean(nsB)   ,std(nsB),  min(nsB), max(nsB) }
//    平均    標準偏差   最小   最大
 ans = {4.99495, 0.30041 , 3.82579, 6.27571}

//理論値 5      0.3

 

において、

histplot([3.5:0.1:6.5],nsB,5)

nsC = nsA + nsB;

{ mean(nsC) ,std(nsC), min(nsC), max(nsC) }
//    平均    標準偏差    最小    最大
ans = {7.9946, 0.317144, 6.84426 , 9.25621}

histplot([6.5:0.1:9.5],nsC,5)

 


一般にこの様な場合、平均 Xc は Xc = Xa + Xb単純平均 、標準偏差σc は σc=√(σa^2 + σb^2) の平方和(偏差平方和)と なります。(理由は数理統計の本に書いてあるので読んでください)

確かめてみると、

mean(nsC)=7.9946 ≒ 2.99965 + 4.99495 = 7.9946
st_deviation(nsC)=0.31714 ≒ √(0.100441^2 + 0.30041^2) = 0.316756


だから、ほぼあいますな。

というわけで、
「平均 3cm,標準偏差 0.1cm」+「平均 5cm,標準偏差 0.3cm」= 「平均 8cm 標準偏差 0.32cm」

なんですが、これ工学な実生活ではどの様に扱うかというと、

例えばモノの規格で 3cm±0.3cm とかありますね。
分布が正規分布の場合、数理統計学によりますと平均±3σは全体の99.7%をカバーします
ちなみに、Excelの場合は、NORMSDIST()という関数があり、両側値の内部の割合は
=(1-(1-NORMSDIST(sigma))*2)*100 [%]
で求まります。この計算によると

±3σ 99.73%    不良率は 1千個に3個(=昔から、千に三つしか合わないといいますが、まさに言いえて妙でしょう)
±4σ 99.9937%  不良率は 2万個に1個
±5σ 99.999943% 不良率は200万個に1個

なので、まぁ、例えば±3σで考えると  平均3cm , σ=0.1cm だと 規格 3±0.3cm は99.7%カバーするはずです。
それでは nsA の中で3±0.3に入っている割合を数えてみましょう。

fA=find ( (3-0.3) <= nsA & nsA <= (3+0.3) ); //findで 3±0.3 の中に入っている要素のリスト を得る

length(fA) / length(nsA) * 100.0 //全体の数で割る。

  ans = 99.69

    となり、ほぼ計算どおりですね

ちなみに ±4σだと

fA=( (3-0.4) <= nsA & nsA <= (3+0.4) );
length(fA) / length(nsA) * 100.0
  ans = 100

    となって、1万個全部含まれてしまいました。(これも計算どおり)


(従って、現在 工程能力Cp ≡ 規格幅 / 3σ → 4σ/3σ= 1.33  というのが、殆ど不良を出さない工程能力値だとされています。一応

さて、例えばこの3cm±0.3cmっていうのは ±3σ=±0.3cmとしましょう。
すると
「平均 3cm,標準偏差 0.1cm」+「平均 5cm,標準偏差 0.3cm」= 「平均 8cm 標準偏差 0.32cm」
っていうのは
「 3cm±0.3cm」+「 5cm±0.9cm」= 「 8cm ± 0.96cm」
となります。
単純に算術平均すると
「 3cm±0.3cm」+「 5cm±0.9cm」= 「 8cm ± 1.20cm」
となるんですがね。

さっきの例で、nsCが 8cm ± 0.96cm に入っているかどうか確かめてみると、

fC=( (8-0.96).<= nsC & nsC <= (8+0.96) );
length(fC) / length(nsC) * 100.0
  ans = 99.78

となりますな。

だから実際には
「 3cm±0.3cm」+「 5cm±0.9cm」= 「 8cm ± 1.20cm」までバラつくことは殆どなくて
「 3cm±0.3cm」+「 5cm±0.9cm」= 「 8cm ± 0.96cm」くらいになる ってことですね。

逆に考えれば、±0.96cmに押さえるための材料は±0.3cmと±0.9cmで、ほぼOKって事ですね。
一般に材料の精度を上げれば価格が高くなりますから、このへんのバランスを考えるのが 実務的でしょう。


★正規分布の掛け算

じゃあ、合計じゃなくて、結果が乗算になるときは、どうだろ?例えば縦横の長さが正規分布した時、面積の分布は?

単純に標準正規分布(平均が0、分散が1)で考えてみる。

SCILABのSCIスクリプト : normal.sci はこんな感じで、ちょっとデモっぽく書いてみる。

>>exec( "noemal.sci")

N(0,1) Y mean= 0.01 std= 1.00
N(0,1) Z mean= 0.01 std= 1.00

こんな感じの分布になっている。(赤のNormalっていう線は正規分布関数 :


 

正規分布を足してみると ***

N(0,1)X=Y+Z mean= 0.02 std= 1.41

やっぱり正規分布です。標準偏差が√(1*1+1*1)=√2=1.414なのでちょっと裾野が広がった感じ。

じゃぁ、掛け算はどうだろう

正規分布の掛け算
:N(0,1)X=YxZ mean=-0.01 std= 1.01

というわけで、平均=0、標準偏差=1なんだけど、中央が高くなってしまってますな。

というわけで、 正規分布の掛け算は正規分布にはならない様です。

return