Scilabで数値微分 その1

数値計算から関数の微分を求めようと考えたとき、微分の定義の通りに差分を計算すれば良いということはすぐに思いつきます。

f'(x) = \lim_{h \to 0} \frac{f(x+h)-f(x)}{h}

上記の式のhを次第に小さくしていけばf'(x)の近似値が得られます。
しかしながらhの数を小さくしすぎると数値計算の丸め誤差の影響で逆にf'(x)の値からずれていってしまうことが知られています。

これに関連する重要なポイントをまとめると以下の3つになります。

  • 倍精度浮動小数点数を使い丸め誤差を小さくする
  • 打切り誤差の小さいアルゴリズムを使う
  • 丸め誤差と打切り誤差のトータルで誤差が小さくなるhを探す


Scilabは特に指定をしなくても倍精度浮動小数点数演算をしてくれるので、今回は数値計算の常識に倣って数値微分の式とhの大きさを色々と変更したときに誤差の大きさがどのようになるのかを見てみます。


数値微分の近似式


数値微分を計算する公式として数値計算の常識に倣って前進差分f1(x)、中心差分f2(x)、及びそれらのRomberg 1段公式を用いました。

前進差分
f_1(x,h) = \frac{f(x+h)-f(x)}{h}

中心差分
f_2(x,h) = \frac{f(x+h)-f(x-h)}{2h}

前進差分に対するRomberg 1段
\frac{f_1(x,h)-\frac{1}{2}f_1(x,2h)}{1-\frac{1}{2}}

中心差分に対するRomberg 1段
\frac{f_2(x,h)-\frac{1}{4}f_2(x,2h)}{1-\frac{1}{4}}

Scilabプログラム


実際に計算する関数は普通のサイン関数f(x)=sin(x)とします。
当然ながら解析的に微分が可能でf'(x)=cos(x)です。

数値計算の常識と同様にx=0.3πのときの微分を計算します。

微分の刻み幅hは
h = 2^{-n}

とし、nが大きくなるほどhは小さくなります。

作成したScilabのプログラムはdiff_sce.txtです。

clear;

// *** 刻み幅の設定 ***
n = [0:1:50];
h = 2 ^ (- n);

x = 0.3 * %pi;

// *** 微分の近似値 ***
// 前進差分
f1 = (sin(x + h) - sin(x)) ./ h;
// 中心差分
f2 = (sin(x + h) - sin(x - h)) ./ (2 * h);
// 前進差分に対するRomberg 1段公式
romberg1 = 2 * ((sin(x + h) - sin(x)) ./ h - 0.5 * (sin(x + 2 * h) - sin(x)) ./ (2 * h));
// 中心差分に対するRomberg 1段
romberg2 = ((sin(x + h) - sin(x - h)) ./ (2 * h) - 0.25 * (sin(x + 2 * h) - sin(x - 2 * h)) ./ (4 * h)) / 0.75;

// *** グラフの軸設定 ***
a = gca();
a.data_bounds = [0,1E-16;50,1];
a.log_flags = "nl";

// *** グラフのプロット ***
plot(n,abs(f1 - cos(x)),'-sr');
plot(n,abs(f2 - cos(x)),'-sm');
plot(n,abs(romberg1 - cos(x)),'-sb');
plot(n,abs(romberg2 - cos(x)),'-sg');

// *** グラフの体裁 ***
xlabel("n");
ylabel("Err");
legend(['f1(x)','f2(x)','Romberg f1(x)','Romberg f2(x)'],4);


結果


以下に示すFig.1は上記の4つの式で計算したf(x)のx=0.3πでの微分の値とf'(x)=cos(x)の値の差の絶対値を取ったものです。
この2つの差が数値計算の誤差ということになります。

001_20131110042107c73.png

Fig.1: 数値微分と解析的な微分の間の関係。縦軸は誤差の絶対値、横軸はnで、nが大きくなるほどhは小さくなる(h = 2-n)。


まずn>25では、全てのグラフにおいて同じ挙動を示しています。これが丸め誤差によるものでどの式を使っても丸め誤差の問題は改善されないことが分かります。

次にnが小さい側、すなわち打切り誤差の影響を見てみます。
前進差分f1(x)(赤)と中心差分f2(x)(ピンク)を比べてみると、nが小さいとき中心差分f2(x)(ピンク)の方が誤差が小さくなることが分かります。

また前進差分f1(x)(赤)とそのRomberg 1段(青)を比較すると、Romberg 1段(青)で打切り誤差が改善していることが見て取れます。

中心差分f2(x)(ピンク)とf1(x)に対するRomberg 1段(青)では誤差はほとんど変わりませんが、中心差分f2(x)(ピンク)にさらにRomberg 1段(緑)をおこなうとこれら4つの式の中で最も誤差の少ない式に出来ることがわかります。

関連エントリ




付録


このエントリで使用したScilabのプログラムを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 数値微分 

緑色のマグカップ

先日twitterにも投稿したのですが、近所の東急ストアでマグカップを購入しました。( https://twitter.com/gomisai/status/446610830826536960 )

DSCF1641.jpg


実用性に関しては、なんとも言えない感じですが、こういうデザインは私は大好きです(笑)。
緑色をしているので緑茶を入れると様になります。紅茶やウーロン茶もありです。でも、カルピスを入れると微妙な気分に。

Scilabで数値積分:地球深部の密度と圧力

Scilabで数値積分: 固体の比熱ではintegrate関数を用いてデバイモデルから結晶固体の比熱を計算しました。
またScilabで乱数による関数の積分ではモンテカルロ法を用いて関数の積分を行いました。

これらはいずれも既知の関数f(x)について具体的な積分範囲、たとえばa≦x≦bでの積分値

\int_a^b f(x)\mathrm{d}x

を求める手法でした。

これらの方法では関数f(x)の形が式として分かっている必要があります。
しかしながら、離散的な数値データに対する積分値が計算したい場合もあります。

001_201311020335151bb.png

Fig.1: 地球の内部構造(世界初!地球中心部の超高圧高温状態を実現 ~ようやく手が届いた地球コア~より)


今回は地球内部の密度を表す数値表から、数値積分を使って重力加速度と圧力の深さ依存性を計算します。



標準地球モデル:PREM


地球内部の構造は、地震波観測からfig.1のような層構造をしていることが知られています。
もっとも有名な地球内部構造のモデルはDziewonski and Anderson (1981)によるPREM(Preliminary Reference Earth Model)と呼ばれるものです。

002_20131102033515160.png

Fig.2: 地球内部の地震波速度分布と密度分布。横軸は地球の中心からの距離(半径)。


PREMにはP波速度、S波速度、密度の深さごとの値、及びこれら3つの数値から計算できる重力加速度、圧力などが表として掲載されています。(参考:PREMの表地球を貫通するトンネル内の落下運動)

今回はPREMの密度から重力加速度と圧力を計算し、PREMの論文の中の値と比較して計算結果のチェックを行います。

重力加速度の計算


高校物理では重力加速度は約9.8[m/s]として計算をすると思います。
重力加速度はは大きな質量を持つ地球から万有引力によって引っ張られる力に起因しているため、地球内部では重力加速度の値が地球表面の値とは異なります。

地球の中心からr[m]での密度をρ(r)[g/m3]とするとr[m]より深い部分の質量M(r)[g]は

M(r) = 4 \pi \int_0^r u^2 \rho(u) \mathrm{d}u

重力加速度は、万有引力の法則より

g(r) = G\frac{M(r)}{r} = \frac{4\pi G}{r}\int_{0}^{r}u^2\rho(u)\mathrm{d}u

ここでGは万有引力定数G=6.67259×10-14[m3/s2/g]です。
PREMの数値データはPREM.txtに保存しfscanfMat等を用いてScilabから読み出しを行います。離散データの数値積分にはinttrapを用いました。

重力加速度計算の積分のところだけ抜き出すと、以下のようになります。

// *** 重力加速度の計算 ***
Gravity_num = zeros(Radius);
for i = 2:length(Gravity_num)
Gravity_num(i) = 4 * %pi * grav / (Radius(i) ^ 2) * inttrap(Radius(1:i), Radius(1:i) .^ 2 .* RHO(1:i));
end


圧力の計算


地球内部では非常に長い地質学的時間のなかで静水圧平衡が成り立っていると考えることが出来るため、圧力の計算が出来ます。

\frac{\mathrm{d}P}{\mathrm{d}r} = - \rho(r)g(r)

これを素直に積分すると以下の式が得られます。

\int\mathrm{d}P = - \int\rho(r)g(r)\mathrm{d}r

文字を置き換えて

\int\mathrm{d}Q = - \int\rho(u)g(u)\mathrm{d}u

地球表面rsでの圧力Ps(つまり大気圧)から地球深部rの圧力Pまで積分すると

\int_{P_s}^{P}\mathrm{d}Q = - \int_{r_s}^r\rho(u)g(u)\mathrm{d}u

大気圧はPs=1013.25[hPa]ですが、これは地球内部の圧力である数[GPa](数万気圧)と比べて充分小さいので左辺はほぼ地球深部の圧力Pと考えることが出来ます。

\int_{P_s}^{P}\mathrm{d}Q = P - P_s \simeq P

したがって圧力は以下のような積分で計算することが出来ます。

P = \int_{r}^{r_s}\rho(u)g(u)\mathrm{d}u

この積分を行う部分をScilabのプログラムから抜粋すると以下のようになります。

// *** 圧力の計算 ***
Pressure_num = zeros(Radius);
for i = 1:length(Pressure_num)-1
Pressure_num(i) = inttrap(Radius(i:$), RHO(i:$) .* Gravity_num(i:$));
end


計算結果


重力加速度と圧力の計算結果をFig.3に示します。

003_20131102033515738.png

Fig.3: 地球内部の重力加速度(上)と圧力(下)


赤の実線で示されたのがPREMの表にある値で、青の破線で示したのが今回積分から計算された値です。多少の誤差が発生しているようですが、おおよそ上手く計算できています。

全ソースコードはDeepEarth_sce.txtです。

clear;

// *** 物理定数 ***
grav = 6.67259e-14; // 万有引力定数 (m^3/s^2 g)

// *** PREMのテーブルを読み出し ***
X = fscanfMat('PREM.txt');
Radius = 1E3 * X(:,1); // 半径 (m)
Vp = X(:,2); // P波速度 (m/s)
Vs = X(:,3); // S波速度 (m/s)
RHO = 1E3 * X(:,4); // 密度 (g/m^3)
Ks = 1E12 * X(:,5); // 断熱体積弾性率(Pa)
Mu = 1E12 * X(:,6); // 剛性率(Pa)
Nu = 1E12 * X(:,7); //
Pressure = 1E12 * X(:,8); // 圧力(Pa)
Gravity = X(:,9); // 重力加速度 (m/s^2)

// *** 重力加速度の計算 ***
Gravity_num = zeros(Radius);
for i = 2:length(Gravity_num)
Gravity_num(i) = 4 * %pi * grav / (Radius(i) ^ 2) * inttrap(Radius(1:i), Radius(1:i) .^ 2 .* RHO(1:i));
end

// *** 圧力の計算 ***
Pressure_num = zeros(Radius);
for i = 1:length(Pressure_num)-1
Pressure_num(i) = inttrap(Radius(i:$), RHO(i:$) .* Gravity_num(i:$));
end

// *** 地震波速度と密度のプロット ***
scf(0)
// 地震波速度
xsetech([0,0,1,0.45]);
plot(1E-3 * Radius, 1E-3 * Vp,'-r');
plot(1E-3 * Radius, 1E-3 * Vs,'-b');
legend(["Vp","Vs"],1);
ylabel("Velocity (km/s)");
xlabel("Radius (km)");
// 密度
xsetech([0,0.45,1,0.45]);
plot(1E-3 * Radius, 1E-3 * RHO,'-r');
ylabel("Density (kg/s^3)");
xlabel("Radius (km)");

// *** 重力加速度と圧力のプロット ***
scf(1)
// 重力加速度
xsetech([0,0,1,0.45]);
plot(1E-3 * Radius,Gravity,'-r');
plot(1E-3 * Radius,Gravity_num,'--b');
legend(["Table","Numerical"],4);
ylabel("Acceleration of gravity (m/s^2)");
xlabel("Radius (km)");
// 圧力
xsetech([0,0.45,1,0.45]);
plot(1E-3 * Radius,1E-12 * Pressure,'-r');
plot(1E-3 * Radius,1E-12 * Pressure_num,'--b');
legend(["Table","Numerical"],1);
ylabel("Pressure (GPa)");
xlabel("Radius (km)");


台形積分とスプライン積分


今回使用したinttrapは、データ店の間を直線で補完する(要するに台形の面積を計算する)台形積分です。
もう少し高度な計算をするスプライン積分もintsplinを使って計算することが出来ます。
入力するパラメータも同じなので、多くの場合そのまま置き換えることが出来るでしょう。

ただし今回のPREMのようにデータが不連続となる点で複数の値を持っているような表の積分を行うときにはintsplinでは下記のようなエラーが出ます。

!--error 999 
splin: 入力引数 #1 の値が間違っています: (厳密に)昇順ではありません, もしくは +-inf を検出しました.


このような場合は、今回のように台形積分inttrapを使います。

関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 数値積分  マントル 

Scilabで円周率の近似値

3月14日は円周率の日だからというわけではないですが、Scilabで楽しむ確率論に倣ってモンテカルロ法を用いた円周率の近似値計算を行いました。

001_20131105060858590.jpg

Fig.1: πのパイ(Wikipedia Pi pie2.jpg)


半径1の円の面積はπなので、その四分の一の面積はπ/4です。
したがって0≦x≦1, 0≦y≦1の範囲にランダムに点を打ち、半径r=1の円の中に入っている数を数えれば、円周率の近似値を求めることが出来ます。

今回はn=10000個の乱数の組(x,y)が半径r=1の円の中に入る数から円周率を求め、更にこれをm=1000回繰り返したときどの程度のばらつきがあるのかをヒストグラムにしました。


モンテカルロ法による円周率計算


0≦x≦1, 0≦y≦1の領域に一様分布に従うように点を打ったとき、半径r=1の円の中に入る点の数がπ/4に比例します。
以下のようなScilabプログラムにて円周率の近似値を得ることが出来ます。

n = 10000;
x = rand(1,n);
y = rand(1,n);

// グラフの縦横比を等しくする
h = scf(0); // ウィンドウを作成
clf(); // 表示をクリア
ha = h.children(1); // Axes(座標軸)オブジェクトへのハンドルを取得
ha.isoview = "on"; // 座標軸の縦横比を等しくする
ha.data_bounds = [0, 0; 1, 1]; // 座標軸表示範囲の設定

// 乱数の組(x,y)をプロット
plot(x,y,'.r','markersize',1);
// しきい値となる半径1の円のプロット
xx = linspace(0,1);
plot(xx,sqrt(1 - xx .^ 2),'-b');
// グラフのラベル
xlabel("x");
ylabel("y");

// 1回の計算から得られた円周率
p1 = sum(bool2s(x ^ 2 + y ^ 2 <= 1.0)) / n * 4


002_20131105060857b76.png

Fig.2: [0,1]×[0,1]の領域内にランダムに点を打ったとき、半径1の領域に打たれる点の頻度が近似的にπ/4になる。


繰り返し計算


以上の計算を更にm=1000回繰り返してどの程度のばらつきがあるのかをヒストグラムにしました。
乱数を格納するために大きなメモリを確保する必要があります。これにはstacksize関数を利用します。

m=1のときの計算も含めたScilabプログラムはcircle_sce.txtとなりました。

003_2013110506085769a.png

Fig.3: (n,m)=(10000,1000)のときのヒストグラム


関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab モンテカルロ解析 乱数 ヒストグラム 円周率 

Scilabで乱数による関数の積分

Scilabには関数の積分の計算を行うintegrateがあるため簡単な関数の積分には必要が無いのですが、乱数を使った積分の方法について書きます。
[0,1]の範囲で一様分布な乱数列を{Xk}k=1,2...nとすると関数f(x)の積分は以下のようにして計算することが出来ます。

\begin{eqnarray}\int_a^b f(x) \mathrm{d}x & = & \int_0^1 f(a+(b-a)r)(b-a)\mathrm{d}r \\ & = & \lim_{n \to \infty} \frac{1}{n}\sum_{k=1}^{n}f(a+(b-a)X_k)(b-a) \end{eqnarray}


乱数による積分の考え方


[0,1]の範囲の一様分布な乱数列{Xk}k=1,2...nはScilabを使うと簡単に得られます。(参考:Scilabで乱数の生成,Scilabでコイン投げ)

このとき関数f(x)の積分は、f(Xk)の平均から求めることが出来るでしょう。

\int_{0}^{1}f(x)\mathrm{d}x = \lim_{n \to \infty} \frac{1}{n} \sum_{k=1}^{n}f(X_k)

例として、以下に挙げる具体的な関数f(x)に対して積分の計算をします。

f(x) = \exp(-x^2)

関数f(x)のグラフはFig.1のような形になります。積分は0≦x≦1の範囲の面積を求めることと同じです。


001_20130926022100898.png
Fig.1: 関数f(x)=exp(-x2)のグラフ


Scilabで楽しむ確率論(PDF)では、本文中に書いてある数式と実際にScilabで計算している数式が異なっています。本エントリではプログラムの方にあわせることにします。

モンテカルロ法以外の計算法


実際のところ、この程度の関数なら乱数を用いたモンテカルロ法を使うよりも、素直に数値積分をするほうがはるかに良い結果が得られます。

Scilabではintegrateをつかって簡単に数値積分が出来ます。(参考:Scilabで数値積分: 固体の比熱)
また、もっと基本的な方法として、解析的に原始関数を決められないか考えるという方法もあります。Maximaを使って不定積分を行うと以下のような結果が得られました。

\int f(x)\mathrm{d}x = \frac{\sqrt{\pi} \mathrm{erf}(x)}{2} + C(C: 積分定数)

まあこの場合は誤差関数erfが入っているので数値積分から逃れられてはいないというのも事実ではありますが。

数値積分、解析的な積分、乱数を使った積分のそれぞれの方法で積分値を計算するScilabのプログラムがfunc_sce.txtです。
数値積分と解析的な積分の両方の方法で、積分の値は0.7468241となりました。モンテカルロ法で求めた値もおおよそ同じくらいになります。

clear;

// *** 関数の定義 ***
// 関数f(x)
deff('y = f(x)','y = exp(- x .^ 2)');
// 原始関数F(x)
deff('y = F(x)','y = sqrt(%pi) .* erf(x) ./ 2');

X = linspace(0,5);

// *** グラフの描画 ***
plot(X,f(X));
xlabel("x");
ylabel("f(x)");

// *** 積分の計算 ***
// 解析的な積分
F(1) - F(0)
// 数値積分
integrate('f(x)','x',0,1)
// モンテカルロ積分
n = 1000;
mean(f(rand(1,n)))


ヒストグラム


このモンテカルロ法による積分を複数回行ったときのヒストグラムをFig.2へ示します。
確かにintegrateで求めた値を中心にばらついていることが分かります。


002_20130926022059bb7.png
Fig.2: 積分結果のばらつきを示すヒストグラム


積分範囲の変更


次に積分範囲がa≦x≦bの場合へ拡張します。方法は単純モンテカルロ積分 (2)の通りです。

a \leq x \leq b
0 \leq x - a \leq b - a
0 \leq \frac{x - a}{b - a} \leq 1

\frac{x-a}{b-a} = r

とおくと0≦r≦1となるので、置換積分を考えればよいことになります。

x = a+(b-a)r
\frac{\mathrm{d}x}{\mathrm{d}r} = b-a
\mathrm{d}x=(b-a)\mathrm{d}r

したがって

\begin{eqnarray}\int_a^b f(x) \mathrm{d}x & = & \int_0^1 f(a+(b-a)r)(b-a)\mathrm{d}r \\ & = & \lim_{n \to \infty} \frac{1}{n}\sum_{k=1}^{n}f(a+(b-a)X_k)(b-a) \end{eqnarray}

となります。

関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 数値積分 乱数 モンテカルロ解析 Maxima 

Scilabでコイン投げ

Scilabで楽しむ確率論(PDF)に倣ってコイントスのシミュレーションをScilabを用いて行いました。
その結果、コイントスの確率変数が0.5に収束することが確認できました。また、このときの分布が正規分布になることも確認できました。


001_2013092520034186d.png
Fig.1: 確率変数の収束



確率変数の収束


Scilabで楽しむ確率論(PDF)のコイン投げのシミュレーションを行います。Scilabで乱数の生成で調べたとおりrandは0~1の間の一様乱数を発生させます。

コイントスをしたとき、必ず表か裏が2分の1の確率で出るので、生成された乱数の値が0.5より小さいとき表が出たと考え、表と裏をそれぞれ数値の1と0に対応させると考えます。
k回目のコイントスの結果をXkとすると1回目からn回目の結果を並べた確率変数列{Xk}k=1,2,...nを考えることが出来ます。

ここでŜnを以下のように定義します。

\hat{S}_n = \frac{1}{n}\sum_{k=1}^n X_k

これはn回コインを投げて表が出る「頻度」を表しています。
例えば3回コインを投げて表(1)、裏(0)、裏(0)と出たとします。

\hat{S}_3 = \frac{1}{3}(1 + 0 + 0) = 0.333...

表が出るか裏が出るかはランダムなのでŜnの値もコイントスの結果によって変わります。しかしながらnが大きくなるに従ってŜnが0.5に近づくであろう事は直感的に予想できます。
これをScilabを用いて1000回(n=1000)のコイントスについて実際に計算したものがFig.1です。

実際にn=1000のときŜnが0.5に近くなっていることが確認できます。
このことを確率変数の収束と呼びます。

中心極限定理


次に 『1000回コイントスを行う(n=1000)』 という事を5回繰り返す(m=5)ということを考えます。
その結果がFig.2です。


002_20130925200341118.png
Fig.2: Ŝnのばらつき


それぞれの線のn=1000のときの値は予想通り0.5に近くなってはいますが、全ての結果が完全に0.5になっているわけではなく、値にばらつきがあります。

『n回のコイントス』をm回繰り返したときのŜnのばらつきは正規分布に従うことが知られており、このことを中心極限定理と呼びます。

標準正規分布関数

f(x) = \frac{1}{\sqrt{2 \pi}}\exp \left(-\frac{x^2}{2}\right)

と比較するためにŜnを適切にシフト・スケーリングするとZnが得られます。

Z_n = \frac{n(\hat{S}_n - p)}{\sqrt{np(1-p)}}

n=1000, m=1200 についてZnのヒストグラムを標準正規分布関数と比較したものがFig.3です。


003_20130925200340892.png
Fig.3: Znのヒストグラムと標準正規分布関数


Appendix: 正規分布


LTspiceモンテカルロ解析の定数分布 その2で抵抗器の部品定数の分布は正規分布に従うと書き、秋月電子通商で購入したカーボン抵抗の平均μや分散σ2を求めました。

正規分布の確率密度関数は以下のように書くことが出来ます。

f(x) = \frac{1}{\sqrt{2 \pi \sigma^2}}\exp \left(-\frac{(x - \mu)^2}{2 \sigma^2}\right)

この正規分布をN(μ,σ2)と書きます。
μ=0,σ2=1のときN(0,1)は特別に標準正規分布と呼ばれます。

f(x) = \frac{1}{\sqrt{2 \pi}}\exp \left(-\frac{x^2}{2}\right)

Scilabで楽しむ確率論(PDF)では「n回のコイントスをN回繰り返す」となっていますが、正規分布N(μ,σ2)と紛らわしいので本エントリでは「n回のコイントスをm回繰り返す」としました。

関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab モンテカルロ解析 確率論 乱数 

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoCOPアンプCPA強磁性PICモンテカルロ解析常微分方程式odeトランジスタecalj状態密度DOSインターフェース定電流スイッチング回路PDS5022半導体シェルスクリプト乱数レベルシフトHP6632A温度解析分散関係I2Cトランジスタ技術R6452A可変抵抗ブレッドボードセミナーバンドギャップ数値積分確率論反強磁性偏微分方程式バンド構造絶縁熱設計非線形方程式ソルバフォトカプラシュミットトリガLEDLM358カオスISO-I2C三端子レギュレータGW近似A/Dコンバータカレントミラーアナログスイッチ数値微分マフィンティン半径TL431発振回路サーボPC817CUSB直流動作点解析74HC4053補間FFTBSch開発環境パラメトリック解析2ちゃんねるチョッパアンプ量子力学bzqlty電子負荷イジング模型LDA標準ロジックアセンブラ基本並進ベクトルブラべ格子単振り子熱伝導位相図TLP621キュリー温度繰り返し状態方程式MaximaVESTAスイッチト・キャパシタ相対論FETランダムウォークスピン軌道相互作用SMP六方最密充填構造抵抗不規則合金ewidthスレーターポーリング曲線GGAラプラス方程式cygwingfortranQSGW失敗談コバルト条件分岐TLP521テスタLM555Writer509TLP552格子比熱マントルデータロガー自動計測詰め回路ガイガー管ダイヤモンドQNAPMCUFXA-7020ZR過渡解析三角波UPSNE555固有値問題熱力学ブラウン運動フェルミ面awk起電力第一原理計算OpenMPfsolveubuntu最大値xcrysden最小値最適化仮想結晶近似VCA差し込みグラフスーパーセル井戸型ポテンシャル平均場近似シュレディンガー方程式FSMフラクタルOPA2277固定スピンモーメント2SC1815全エネルギー合金multiplotgnuplotc/aTeX結晶磁気異方性interp1ウィグナーザイツ胞初期値マンデルブロ集合疎行列面心立方構造fcc不純物問題非線型方程式ソルバフィルタL10構造PGA半金属二相共存SICZnOウルツ鉱構造BaO重積分クーロン散乱磁気モーメント電荷密度三次元CIF岩塩構造CapSenseノコギリ波デバイ模型ハーフメタル正規分布フォノンquantumESPRESSOルチル構造スワップ領域リジッドバンド模型edelt縮退キーボード軸ラベルグラフの分割凡例トラックボールPC不規則局所モーメント片対数グラフトランス両対数グラフCK1026MAS830L直流解析Excel円周率パラメータ・モデルヒストグラム日本語最小二乗法等価回路モデルGimp線種シンボルTS-110TS-112PIC16F785LMC662化学反応文字列specx.f入出力ifortマテリアルデザインヒストグラム確率論Realforce等高線ジバニャン方程式P-10Ubuntuナイキスト線図Crank-Nicolson法陰解法熱拡散方程式HiLAPWAACircuit連立一次方程式負帰還安定性境界条件EAGLEMBE関数フィッティング

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ