←こんなやつで、ＳＣＩＬＡＢの正規分布乱数を使ったモンテカルロシミュレーション
作成：2000-04-06
更新：2000-10-03
★正規分布の足し算
正規分布っていうのは←こんなやつで、
←こんなやつですが、
普通の乱数っていうのは０〜１まで均等に分布しますが、工学のシミュレーションをするときには、だいたい、モノの精度って正規分布してますから、正規分布な乱数が欲しいですね。
技術評論社：奥村晴彦著：C言語によるアルゴリズム辞典 によれば rnd()を0〜1.00の一様乱数として、
t=sqrt(-2.0 * log(1-rnd()));
u=2*PI*rnd();
r1=t*cos(u) ; r2=t*sin(u);
で２個ずつ、平均０、分散１の正規分布する乱数が得られるそうです。
SCILABも（多分）これによる正規分布乱数を発生する関数 rand(,,"normal")っていうのが実装されてますので早速使ってみましょう。
//平均０、標準偏差１の正規分布乱数は randn() で求められる。 この乱数を10000個発生させて、
|
ns = rand(1,10000 , "normal")； //
1×10000 個の行列＝１次元配列 mean(ns) //平均 st_deviation(ns) 標準偏差 min(ns) 最小値 max(ns) 最大値 |
ヒストグラムを書いてみよう。
って簡単にヒストグラムが書けるところが便利っすね
|
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 というＭａＴＸとの記述とは違う事に注意（この辺、ツールで異なるのはしかたがないんだが、なんとかならんもんか？）
この場合、
|
mean(ns) //平均 st_deviation(ns) //標準偏差 min(ns) //最小値 max(ns) //最大値 |
んじゃ、応用問題で、
長さ 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; { mean(nsA) ,st_deviation(nsA), min(nsA), max(nsA) } { mean(nsB) ,st_deviation(nsB),min(nsB), max(nsB) } //この nｓＡ と nｓＢ から nｓＣ を作る。 nsC = nsA + nsB; { mean(nsC) ,std(nsC), min(nsC), max(nsC) } |
|
{ mean(nsA) ,std(nsA), min(nsA), max(nsA) } //理論値 3 0.1
|
histplot([1.5:0.1:4.5],nsA,5) |
|
{ mean(nsB) ,std(nsB), min(nsB), max(nsB) } //理論値 5 0.3
において、 |
histplot([3.5:0.1:6.5],nsB,5) |
|
nsC = nsA + nsB; { mean(nsC) ,std(nsC), min(nsC), max(nsC) } |
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 とかありますね。
分布が正規分布の場合、数理統計学によりますと平均±３σは全体の99.7%をカバーします
ちなみに、Excelの場合は、NORMSDIST()という関数があり、両側値の内部の割合は
=(1-(1-NORMSDIST(sigma))*2)*100 [%]
で求まります。この計算によると
±3σ 99.73% 不良率は 1千個に３個（＝昔から、千に三つしか合わないといいますが、まさに言いえて妙でしょう）
±4σ 99.9937% 不良率は 2万個に１個
±5σ 99.999943% 不良率は200万個に１個
なので、まぁ、例えば±３σで考えると 平均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) ); |
となって、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) ); |
となりますな。
だから実際には
「 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って事ですね。
一般に材料の精度を上げれば価格が高くなりますから、このへんのバランスを考えるのが 実務的でしょう。
じゃあ、合計じゃなくて、結果が乗算になるときは、どうだろ？例えば縦横の長さが正規分布した時、面積の分布は？
単純に標準正規分布（平均が０、分散が１）で考えてみる。
SCILABのSCIスクリプト : normal.sci はこんな感じで、ちょっとデモっぽく書いてみる。
|
>>exec( "noemal.sci") N(0,1) Y 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なのでちょっと裾野が広がった感じ。 じゃぁ、掛け算はどうだろう？ 正規分布の掛け算
というわけで、平均＝０、標準偏差＝１なんだけど、中央が高くなってしまってますな。 というわけで、 正規分布の掛け算は正規分布にはならない様です。 |
