正規分布

正規分布の確率密度関数を描き、範囲を指定して確率を計算します。

よくあるのは、平均m、標準偏差σに対し、m±σ , m±２σ , m±３σ に入る確率です。

次の図はm=0の場合です。σ=1 のとき標準正規分布です。±σに入る確率はおよそ 0.68269 と計算されています。

スクリプトは次の通りです。Drawスロットに書きます。

================================================

σ=C.x-A.x;

m=0;

f(x):=(1/(sqrt(2*pi)*σ))*exp(-(x-m)^2/(2*σ^2));

plot(4*f(#)); // ｙ軸方向に4倍して表示

apply(-10..10,drawtext([#-0.1,-0.2],#,size->12));

apply(1..10,drawtext([-0.3,#*4/10-0.05],#/10,size->12));

drawtext([G.x+0.1,G.y],"m="+m,size->18);

drawtext([C.x+0.1,C.y],"σ="+σ,size->18);

// シンプソン公式による数値積分

Z.y=0; // 点Z,Dがｘ軸から離れないようにする。

D.xy=[Z.x,4*f(Z.x)];

sx=0; //定義域左端

ex=Z.x; //定義域右端

if(sx>ex,ex=0;sx=Z.x);

n=10; // 分点の数

dx=(ex-sx)/n; //区間幅

xn=apply(0..n,sx+#*dx); // sxからexまでをn等分した分点のリスト

yn=apply(xn,f(#));

sum=0;

dx=dx/2;

repeat(n,i,

x0=xn_i;

x1=(xn_i+xn_(i+1))/2;

x2=xn_(i+1);

y0=yn_i;

y1=f(x1);

y2=yn_(i+1);

sum=sum+dx*(y0+4*y1+y2)/3;

);

drawtext([1,1.5],"確率="+format(2*sum,5),size->18);

pt=[[Z.x,0],[-Z.x,0],[-Z.x,4*f(-Z.x)]];

nx=abs(Z.x*40);

repeat(nx,s,

x=-Z.x+2*s*Z.x/nx;

pt=append(pt,[x,4*f(x)]);

);

fillpoly(pt,alpha->0.2,color->[1,0,0]);

==========================================

この図から、必要な要素だけをKETCindyで書き出します。

右上のスライダは、σを変更するためのスライダですが、これは書き出す必要はありません。

確率密度関数と、ピンクの領域を囲む線、それぞれの点でしょう。緑の点はZです。この座標を使って3本の線分を表示します。

σの値と確率も書き出します。±σの範囲を斜線で示すことにしましょう。

=============================================

Defvar("s="+text(σ));

Deffun("g(x)", ["regional(y)","y=4*(1/(sqrt(2*pi)*s))*exp(-x^2/(2*s^2))","y"]);

Plotdata("1","g(x)","x");

Listplot("1",[[-Z.x,4*f(-Z.x)],[-Z.x,0],[Z.x,0],[Z.x,4*f(Z.x)]]);

Listplot("2",[[-0.1,2],[0.1,2]]);

Letter([[0,2],"w4","0.5",[1,1.5],"e","確率＝"+format(2*sum,5),[1,2],"e","σ="+text(σ),[Z.x,0],"s2",text(Z.x)]);

Hatchdata("1",["ii"],[["gr1","s"],["sg1","n"]]);

==============================================

実行すると次のようになります。

今度は、σの値を変えたものを表示しましょう。

σを.小さくすると縦長の釣り鐘形になります。

< 戻る >