Microsoft PowerPoint プレゼンテーション.pptx

(1)

Microsoft Excel ®で学ぶ

モンテカルロ計算による

不確かさ評価

本ファイルに含まれるプログラムはモンテカルロ計算の仕組みの理解を促進する目的で私的に使用することを想定したものです。実用的な目的への使用は許可していません。プログラムを使用したことにより損害が発生した場合、作成者やその所属機関は一切責任を持ちません。実務に使用するためのプログラムの開発に興味のある方は、作成者までご連絡下さい。産業技術総合研究所城野克広（最終改訂 2018年12月5日） https://staff.aist.go.jp/k.shirono/

(2)

1. モンテカルロ

計算で不確かさ

評価（原理）

(3)

3

ある球の密度を知るのに、直径の測定と質量の測定を

する。半径の値rは平均1.0 cm、標準偏差0.1 cmの正

規分布に従って測定されたとする。質量の値mは平均

10.0 g、標準偏差0.1 gの正規分布に従うとする。密度

の値は以下の式で値は2.4 g∙cm

−3

_{となる。標準不確か}

さと95 %の信頼の水準の区間を求めよ。

3 3

4

3

4 _r

m

r

m





















g∙cm

−3

(4)

4

不確かさの伝播則を用いて標準不確かさ評価する。



をrで微分すると、

 

₄

7 .

2 g

cm

4

3

3 







r

m





をmで微分すると、

1

₃

0 .

24 cm

-3

4

3 

r



 感度係数

−7.2 g∙cm

−4

(5)

5

 rによる



の標準不確かさ

• rに対する感度係数 −7.2 g∙cm

−4

• rの標準不確かさ0.1 cm

→ これらを乗じて、0.72 g∙cm

−3

_{と計算される。}

 mによる



の標準不確かさ

• mに対する感度係数 0.24 cm

−3

• mの標準不確かさ0.1 g

→ これらを乗じて、0.024 g∙cm

−3

と計算される。





の標準不確かさu(



)

u



0.72 0.024 g∙cm−3 = 0.72 g∙cm−3

(6)

6





の拡張不確かさU

中心極限定理から正規分布として、包含係数２を用

いると、

U = 2×u(



) = 1.4 g∙cm

−3

 95 %の信頼の水準の区間

(



– U,



+ U) = (2.4 – 1.4, 2.4 + 1.4) g∙cm

−3

= (1.0, 3.8) g∙cm

−3

(7)

• バジェットシートは以下のようになる。

不確かさ要因記号標準不確かさ u(x_i) 感度係数 c_i 測定量の標準不確かさ |c_i|u(x_i) (g∙cm−3₎ 備考半径の測定 _u(r) _0.1 cm −7.2 g∙cm−4 0.72 質量の測定 _u(m) _0.1 g 0.24 cm−3 0.024 合成標準不確かさ uc() 0.72 拡張不確かさ U 1.44

k = 2

7

(8)

• モンテカルロ法でやってみる。

• エクセルで標準偏差０、分散１の正規分布は

セルに「=NORM.INV(RAND(),0,1)」と打ち込む

と得られる。

(9)

• 平均1.0 cm、標準偏差0.1 cmの正規分布を

得るには、「 =NORM.INV(RAND(),1.0,0.1) 」と

打ち込む。

(10)

(11)

• 平均10.0 g、標準偏差0.1 gの正規分布を得る

には、「 =NORM.INV(RAND(),10.0,0.1) 」と打ち

込む。

(12)

(13)

• A1に入っている「半径r」、B1に入っている「質

量m」から、「密度r」を算出するために、C1に

「=3/(4PI())B1/A1^3 」と入力します。

(14)

(15)

• 標準偏差を「=STDEV.S(C1:C100)」で計算しま

す。

(16)

16

• 95 %の信頼の水準を持つ区間の最小値は、

全体の2.5パーセント点である。PERCENTILE関

数により、パーセント点は計算できる。この場

合、「=PERCENTILE(C1:C100,0.025)」とする。

(17)

17

• 95 %の信頼の水準を持つ区間の最大値は、

全体の97.5パーセント点である。PERCENTILE

関数により、パーセント点は計算できる。この

場合、「=PERCENTILE(C1:C100,0.975)」とする。

(18)

18

• このデータから、

標準不確かさは0.83 g∙cm

−3

、

95 %の信頼の水準の区間は(1.4, 4.8) g∙cm

−3

と分かる。

• 感度係数を使って計算した95 %の信頼の水

準の区間は(1.0, 3.8) g∙cm

−3

(19)

19

さらに計算を続けると以下のヒストグラムを得る。

観測回数

2.5 %点と 97.5 %点

ヒストグラム：ある定められた区間にあてはまる値が何回現れたかを示す棒グラフ

u(



) = 0.83 g∙cm

–3

(1.4, 4.6) g∙cm

–3



/ g∙cm

–3

(20)

モンテカルロ計算では

1. 感度係数は計算しない。（→ 感度係数が計算

できないときにモンテカルロ計算は便利）

2. 中心極限定理は使わない。（→ 中心極限定理

が使えないときにモンテカルロ計算は便利）

3. タイプAの不確かさの計算が違う。（→ 後ほど）

(21)

21

■便利な点

・非正規性や非線形性が問題になる場合となら

ない場合の場合分けが必要でない。

■問題点

・ソフトウェアの正しさ（ISO 17025 5.4.7a）の証明

をどのようにするのかが一般的でない。

・モンテカルロ法による計算結果を校正結果とし

て受け取った場合の取り扱いがどうしてよいか

分からない。

(22)

22 http://www.bipm.org/utils/common/documents/jcgm/JCGM_101_2008_E.pdf Evaluation of measurement data – Supplement 1 to the "Guide to the expression of uncertainty in measurement" – Propagation of distributions using a

Monte Carlo

method

JCGM 101

(23)

23

2. VBA

プログラミング

入門

(24)

「開発」タブを追加する①

24 １．「ファイル」タブをクリック２．「オプション」をクリック３．「リボンのユーザー設定」をクリック

(25)

「開発」タブを追加する②

25 ４．「メインタブ」のチェックボックスの中から、「開発」をチェック５．「ＯＫ」ボタンをクリック６．「開発」タブが追加された。

(26)

マクロを使う①

26 １．「開発」タブをクリック２．「マクロ」をクリック３．「マクロ名」に適当な名称を入力し、「作成」ボタンをクリックする

(27)

マクロを使う②

27 ４．マクロ編集画面が立ち上がる。５．Sub …の行とEnd Subの間に「Msgbox(“モンテカルロ法始めました。”)」と入力６．編集画面の右上の最小化ボタンをクリックして、編集画面を最小化し、エクセルを表示する。

(28)

28 ７．「マクロ」をクリック８．作成したマクロを選択し、「実行」ボタンをクリックする。８．メッセージが表示される。

マクロを使う③

９．「ＯＫ」ボタンをクリック

(29)

宣言・型

29

Sub Sengen ()

Dim i As Long

Dim x As Double

i = 1

x = 1#

MsgBox (i)

MsgBox (x)

End Sub

i を整数型として扱うことを宣言する。 x を実数型として扱うことを宣言する。 1#は整数としての1ではなく、実数としての1.000000…のことを意味している。 MsgBox()はメッセージボックスに値を表示させる関数。 iに1という値を記録する。

(30)

配列

30

Sub Hairetsu()

Dim x(1 To 3) As Double

x(1) = 1#

x(2) = 2#

x(3) = 3#

MsgBox (x(1))

MsgBox (x(2))

MsgBox (x(3))

End Sub

x(1)、x(2)、x(3)という一連のデータ（配列）として扱うことを宣言。 x(1)に1.0を代入する。以下同様。

(31)

データの入出力①

31

Sub NyuShutsuRyoku1()

Dim x As Double

x = Range("A1").Value

Range("B1").Value = x

End Sub

セルのA1に入っているデータをxに入力する。 xに入っているデータをセルのB1に出力する。

(32)

データの入出力②

32

Sub NyuShutsuRyoku2()

Dim x As Double

x = Cells(2,1).Value

Cells(2,2).Value = x

End Sub

上から２段目、左から１行目のセル、すなわちA2のセルからデータをxに入力する。 xに入っているデータを上から２段目、左から２行目、すなわちB2のセルに出力する

(33)

四則演算

33

Sub Warizan()

Dim x As Double

Dim y As Double

Dim z As Double

x = 1

y = 2

z = x/y

MsgBox (z)

End Sub

加減乗除にはそれぞれ「+」、「‐」、「*」、「/」が対応する。この場合は、xをyで除する。

(34)

34

Sub KansuCos()

Dim x As Double

Dim y As Double

x = 0

y = cos(x)

MsgBox (y)

End Sub

簡単な関数①

「cos()」はcosを返す関数。簡単な関数はVBAで定義されている。（次スライド参照）

(35)

簡単な関数②

35

x^2, x^3, …

x^(1/2), x^(1/3), …

Abs(x)

Sin(x), Cos(x), Tan(x)

Exp(x)

Log(x)

x

2 _, x

3 _, …

√x,

3 √x, …

|x|

sin(x), cos(x), tan(x)

exp(x)

ln(x)

log

₁₀

(x)はLog(x)/Log(10#)

その他にもたくさんの関数が用意されている。

(36)

ワークシート関数①

36

エクセルで使える関数が、VBAではそのままでは

使えないことがある。例えば、最大値を呼び出す

Max関数は、VBAでは準備されていない。しかし、

「WorksheetFunction.Max()」

とすることで、エクセルを用いたVBAではMax関数

を使うことができる。同じように、使用できる関数

が他にもいくつかある。ただし、ワークシート上で

「.」が入っている関数は、それを「_」に置き換える

必要がある場合もある。（具体例はのちほど）

(37)

ワークシート関数②

37

Sub ChooseMax()

Dim x(1 To 3) As Double

Dim y As Double

x(1) = 1#

x(2) = 2#

x(3) = 3#

y = WorksheetFunction.Max(x)

MsgBox (y)

End Sub

x(1) ~ x(3)の中から最も大きいものを選び、yに記録する。このように配列を入力にできる関数もある。

(38)

繰返し（For文）

38

Sub Kurikaeshi()

Dim i As Long

Dim n As Long

n = 0

For i = 1 To 10

n = n + i

MsgBox (n)

End Sub

「For 変数 = 1 To 10」で、その変数を1 から10まで変えながら、「Next」と書いてある行までの演算を繰り返すことを意味している。10を100や1000に変えれば、100回あるいは1000回になる。変数は演算の中で使ってもよいし、使わなくてもよい。後者の場合、変数は単なる繰り返しのための指標となる。「n = n + i」はもともとのnという値にiという値を足して、更新するという意味。例えばnに1、i に2という値が入っているときに、「n = n + i」とすると、nが1から3に更新される。

(39)

演習

39

100から1000までの整数の和はいくつになる

か？これを計算し、出力する「Goukei」というプ

ログラムを作れ。

(40)

条件分岐（If文）①

40

Sub JoukenBunki()

Dim i As Long

i = 1

If (i = 1) Then

MsgBox ("iは1。")

Else

MsgBox ("iは1でない。")

End If

End Sub

「If(（条件式）) Then （演算①） Else （演算②） End If」の構文で条件分岐を行う。（ここで使われた条件式の「i = 1」は「iが1に等しい」という意味。）条件式が真ならば、ThenとElseの間の演算を行う。条件式が偽のとき、Else とEnd Ifの間の演算を行う。もし、偽のときには何の演算もしないならば、Elseを省略して「If(（条件式）) Then （演算①） End If」としてもよい。

(41)

条件分岐（If文）②

41

x = y

x > y

x >= y

x < y

x <= y

x <> y

xとyは等しい。

xはyより大きい。

xはy以上である。

xはyより小さい。

xはy以下である。

xはyではない。

x、yはDouble型でもLong型でもよい。（先のスライド

のように）yの代わりに数字を直接入力してもよい。

(42)

乱数の呼び出し

42

Sub Ransu()

Dim x As Double

x = Rnd

MsgBox (x)

End Sub

「Rnd」は0から1の範囲の一様分布を呼び出す。

(43)

一様分布①

• 不確かさの計算では、最小値と最大値しか分

かっていない値によく用いられる。

x b a (a + b)/2

(44)

一様分布②

44

Sub ItiyouBunpu()

Dim a As Double

Dim b As Double

Dim x As Double

a = ‐100

b = 100

x = (b‐a)*Rnd+a

MsgBox (x)

End Sub

ここでaは一様分布の最小値。ここでbは一様分布の最大値。「(b‐a)*Rnd」とすることで、0から(b‐a) の範囲での一様分布を呼び出せる。それに最小値のaを足すことで、aから bの範囲での一様分布となる。

(45)

ヒストグラム①

• ヒストグラムを作成するためのFrequency関数

を紹介する。まずはワークシートで使ってみる。

45 １．データをあるA列（下の場合A1からA6）に入れる。２．データ区切りの上限となる数値をB列（下の場合、B1からB3）に小さい順に入れる。３．C列のセルを、B列で使ったセルよりも一つ多く選択する。（下の場合C1 からC4）

(46)

ヒストグラム②

4．「=Frequency(A1:A6,B1:B2)」のように Frequency関数に（データのセル）と（データ区切りのセル）の順に入力し、 Cntrlボタン、Shiftボタン、Enterボタンを同時に押す。 46 5．C1のセルにはB1のセルの値以下のデータ数が出力される。C2にはB1のセルの値より大きく、B2のセルの値以下のデータ数が出力される。以下同様で、一番下のセル（下ではC4）にはB列の最大の区切り値より大きいデータの個数が出力される。

(47)

ヒストグラム④

6．C列のデータを選択し、「挿入」タブを選択、「グラフ」から「２－Ｄ縦棒」の「集合縦棒」を選びます。 47 7．グラフが表示される。（実際にこの機能を使うときに、横軸の値は実際の測定値とは関係のない値になっているので、適切に修正する。修正の仕方はのちほど。）

(48)

ヒストグラム⑤

• 次にVBA中で1000個の0から1を範囲とする一

様分布からの乱数がある場合に、ヒストグラ

ムを描く。

48 x １ 0

(49)

ヒストグラム⑥

49

Sub Histgram()

Dim i As Long

Dim x(1 To 1000) As Double

Dim y(1 To 10) As Double

Dim nCount As Variant

For i = 1 To 1000

x(i) = Rnd

(次スライドに続く)

Variantは任意の型ということ。計算の中で自動的に型が与えられる。 Frequency関数の出力はVariant型で指定しないといけない。

(50)

50

For i = 1 To 10

y(i) = i * 0.1

nCount = WorksheetFunction.Frequency(x, y)

For i = 1 To 10

Cells(i, 1) = y(i) ‐ 0.05

Cells(i, 2) = nCount(i, 1)

End Sub

ヒストグラム⑦

nCount(1,1)からnCount(11,1)まで値が格納されているが、nCount(11,1)はゼロなので、出力しないようにした。ここでは、nCount(1,1)が最小値から y(1)の間の配列x中のデータ数、 nCount(2,1)がy(1)からy(2)までのデータ数、以下同様で、nCount(11,1)が y(10)から最大値までのデータ数を返す。上限値ではなく、範囲中の平均値が出力されるように0.05を引いている。

(51)

51 1．出力されたB列のデータを選択し、「挿入」タブを選択、「グラフ」から「２－Ｄ縦棒」の「集合縦棒」を選ぶ。この時点では横軸がデータと合っていない。

ヒストグラム⑧

• ヒストグラムを描いてみよう。

2．「挿入」タブの「データの選択」を選び、。

(52)

52 3．立ち上がったデータソースの選択から、「横軸（項目）ラベル」の「編集」ボタンをクリック。立ち上がった→データ範囲の選択と書いている左の「↑」ボタンをクリック

ヒストグラム⑨

4．この状態でA1からA10を選択すると、下のように選択したセルが入力される。ここで、入力枠の右の「↓」をクリック。

(53)

53

ヒストグラム⑩

5．「軸ラベル」の「OK」ボタンをクリック。続いて、「データソースの選択」の「OK」をクリック。 6．横軸が調整されたヒストグラムを得ることができる。

(54)

正規分布①

54

• 不確かさの計算では、校正結果などの、中心

極限定理がよくあてはまる値に用いられる。



（平均）



（標準偏差）



2 _は分散

(55)

正規分布②

55

• 正規分布から乱数を発生するのに「逆関数

法」を使う。



0

1 x

まずRndで0から1 の一様乱数を得る。確率は0から 1の間の数字である。累積確率分布関数はある分布に 従う値xに対して、確率を返す。 逆に0から1の値を確率と見て、そ れに対応するxを導くこともできる。 この方法である分布の乱数を発生させる方法を逆関数法と呼ぶ。

(56)

正規分布③

56

Sub SeikiBunpu()

Dim mu As Double

Dim sigma As Double

Dim x As Double

mu = 100#

sigma = 10#

x = WorksheetFunction.Norm_Inv(Rnd, mu, sigma)

MsgBox (x)

End Sub

ここでmuは正規分布の平均。 Norm.Inv関数はワークシート上で累積分布関数の逆関数となっている。ここでは、WorksheetFunctionとした上で、「.」を「_」に置き換える必要がある。() の中に、0 ~ 1の一様乱数、平均、標準偏差の順に入力。ここでsigmaは正規分布の標準偏差。

(57)

正規分布④

• これで計算は動く。しかし、ときにRndの値が極めて0あ

るいは1に近いときに、

WorksheetFunction.Norm_Inv(Rnd, mu, sigma)

でうまく計算が行われないことがある。

• これを避けるにはRndが極端に0や1に近い値になるの

を防ぐために、Rndの代わりに

(1‐21.E‐16)Rnd+1E‐16

などとすることがある。

• 「1E‐16」は小さい値で計算に支障のない値なら、他の値

でもよい。このテキストでは、コードの簡明さを重視して

このテクニックは使用しない。

57

(58)

演習

58

標準正規分布から乱数を1000個発生し、その

ヒストグラムを描くための数値を出力するプログ

ラム「Seiki_Hist」を作成せよ。

−５から５の範囲を２０個に区切るものとする。

(59)

t分布①

59

• モンテカルロ計算でタイプＡの不確かさ評価

をするときに使う分布。（通常の不確かさ評価

とは使い方が違うので注意。詳しくは後述。）



= 1



= 3



= 5

は自由度と呼ばれる変数。これを決める とt分布の形が決まる。

(60)

t分布③

60

Sub TBunpu()

Dim nu As Long

Dim x As Double

nu = 4

x = WorksheetFunction.T_Inv(Rnd, nu)

MsgBox (x)

End Sub

ここでも、逆関数法を使う。T.Inv関数はワークシート上で累積分布関数の逆関数となっている。()の中に0 ~ 1の一様乱数、自由度の順に入力。 nuは自由度

(61)

三角分布①

61

• 不確かさ評価では、三角分布の情報が事前

に与えられているときに使う。多くの場合には、

同一の一様分布を２つ重ね合わせた不確か

さ要因である。

(62)

三角分布②

62

• 独立に0 ~ 1の範囲の一様分布に従う２つの

値をx、yとすると、z = x ‒ yは‒ 1 ~ 1の範囲の

三角分布である。

0 x 1 0 ｙ 1 z = x ‒ y 1 ‒ 1

(63)

三角分布③

63

Sub SankakuBunpu()

Dim a As Double

Dim b As Double

Dim x As Double

a = ‐100

b = 100

x = (b‐a)/2*(Rnd‐Rnd)+(b+a)/2

MsgBox (x)

End Sub

ここでaは三角分布の最小値。ここでbは三角分布の最大値。「(b‐a)/2*(Rnd‐Rnd)」とすることで、平均0、全幅(b ‒ a)の三角分布となる。それに平均値の(b+a)/2を足す。

(64)

U字分布①

• 不確かさ評価では、sin関数に従って規則的

に時間変動する変数に対して用いる。

64 0 時間 a _x b x a b

(65)

U字分布②

65

Sub UjiBunpu()

Dim a As Double

Dim b As Double

Dim x As Double

a = ‐100

b = 100

p = WorksheetFunciton.Pi()

x = (b‐a)/2sin(2p*Rnd)+(b+a)/2

MsgBox (x)

End Sub

ここでaとbはU字分布の最小値と最大値。 pには = 3.141592…を代入しておく。 pにを代入した上で、「(b‐a)/2*sin(2*p*Rnd)」とすると、平均 0、全幅(b ‒ a)のU字分布となる。それに平均値の(b+a)/2を足す。

(66)

パーセンタイル点①

• 標準正規分布からの1000個のサンプルを導

き、上側2.5パーセント点を計算する。

66



= 0

?



= 1

2.5 %

(67)

パーセンタイル点②

67

Sub SeikiBunpuPercentile()

Dim i As Long

Dim x(1 To 1000) As Double

Dim y As Double

mu = 0#

sigma = 1#

For i = 1 To 1000

x(i) = WorksheetFunction.Norm_Inv(Rnd, mu, sigma)

y = WorksheetFunction.Percentile(x, 0.975)

MsgBox (y)

End Sub

Percentile関数に（配列、確率）を入力することで、下側のパーセント点を求めることができる。上側100a %点は下側 (100 – 100a) %パーセントである。

(68)

演習

68

自由度9のt分布から、乱数を1000個発生し、上

側の2.5 %点を求めるプログラム「T_Percent」を

作成せよ。

(69)

3. モンテカルロ

計算で不確かさ

評価（応用）

(70)

モンテカルロ計算では

1. 感度係数は計算しない。（→ 感度係数が計算

できないときにモンテカルロ計算は便利）

2. 中心極限定理は使わない。（→ 中心極限定理

が使えないときにモンテカルロ計算は便利）

3. タイプAの不確かさの計算が違う。

(71)

71

• JCGM 101では、

ベイズ統計学

という学

問体系に基づいて、不確かさを計算すること

が基本になっている。

• JCGM 100の従来的な不確かさ評価手法とは

タイプＡの不確かさに対する考え方が違う。

• t分布がタイプAの不確かさに現れるので、

ちょっと面倒に思える一方で、「有効自由度」

を計算する必要がなくなる。

(72)

ある測定値が正規分布に従っている。その平

均



についても、標準偏差



についても事前の

情報はないものとする。４回の測定の結果は、

9、12 、11、10であった。平均



の値の事後分

布はどのように与えられるか。



=_____



= ______

?

72

(73)

正規分布の平均の推定：

正規分布から得た平均の事後確率密度は、繰

り返し数ー１の自由度のt 分布の値に以下をか

けて、平均値を足した値の確率密度である。

残差平方和

繰り返し数×

（繰り返し数－１）

フリークエンティスト統計における平均値の標準偏差 73

(74)

10.5 平均値の標準偏差を計算すると0.65となる。

JCGM 101では、ベイズ統計学に基づき、



は

「自由度４のｔ分布」の値に0.65をかけ、平

均値10.5を足した分布を持つと与える。



の事後確率密度

74 詳細は、 https://staff.aist.go.jp/k.shir ono/downloadBayes.html を参照

(75)

タイプＡの不確かさ評価

75

Sub TypeAUncertainty()

Dim i As Long

Dim n1 As Long

Dim rx1(1 To 4) As Double

Dim mrx1 As Double

Dim srx1 As Double

Dim x1(1 To 1000) As Double

n1 = 4

rx1(1) = 9

rx1(2) = 12

rx1(3) = 11

rx1(4) = 10

(次のページへ続く)

rxが実験データの配列。 xにモンテカルロ法で得たサンプルを記録する。

(76)

76

mrx1 = WorksheetFunction.Average(rx1)

srx1 = WorksheetFunction.StDev_S(rx1) / n1 ^ 0.5

For i = 1 To 1000

x1(i) = mrx1 + srx1 *WorksheetFunction.T_Inv(Rnd, n1 ‐ 1)

MsgBox (WorksheetFunction.StDev_S(x1) )

End Sub

自由度4（= n1 ‐ 1）のt分布の値に計算 した平均値の標準偏差の値(s1)をかけて、平均値(m)を足している。平均と標準偏差はここではワークシート関数を用いて計算する。

(77)

77

ある球の密度を知るのに、直径の測定と質量の測定を

する。半径の値rは平均1.0 cm、標準偏差0.1 cmの正

規分布に従って測定されたとする。質量の値mは平均

10.0 g、標準偏差0.1 gの正規分布に従うとする。密度

の値、標準不確かさと95 %の信頼の水準の区間を求

めよ。（密度の値は以下の式で値は2.4 g∙cm

−3

_{。）}

3 3

4

3

4 _r

m

r

m





















g∙cm

–3

(78)

分布の様子を調べるために、ヒストグラムを作成する。

Sub Density1Histgram()

Dim i As Long

Dim mu1 As Double

Dim sigma1 As Double

Dim x1(1 To 1000000) As Double

Dim mu2 As Double

Dim sigma2 As Double

Dim x2(1 To 1000000) As Double

Dim p As Double

Dim rho(1 To 1000000) As Double

Dim mrho As Double

Dim srho As Double

Dim nCount As Variant

Dim y(1 To 29) As Double

（次スライドに続く）

₇₈ x1に半径の値を入れる。mu1、sigma1 はそれぞれ正規分布の平均と標準偏差。 x2に質量の値を入れる。mu2、sigma2 はそれぞれ正規分布の平均と標準偏差。 rhoに密度の値を入れる。mrhoとsrho は結果として得られる分布の平均と標準偏差。 yはヒストグラムの級の上限値。 nCountにヒストグラムのデータを出力する。

(79)

mu1 = 1#

sigma1 = 0.1

For i = 1 To 1000000

x1(i) = WorksheetFunction.Norm_Inv(Rnd, mu1, sigma1)

mu2 = 10#

sigma2 = 0.1

For i = 1 To 1000000

x2(i) = WorksheetFunction.Norm_Inv(Rnd, mu2, sigma2)

p = WorksheetFunction.pi()

For i = 1 To 1000000

rho(i) = (3 / (4 * p)) * x2(i) / x1(i) ^ 3

（次スライドに続く）

79 rhoに密度の値を入れる。 x1に半径の値を入れる。mu1、sigma1 はそれぞれ正規分布の平均と標準偏差。 x2に質量の値を入れる。 pに = 3.14…を記録する。

(80)

mrho = WorksheetFunction.Average(rho)

srho = WorksheetFunction.StDev_S(rho)

For i = 1 To 29

y(i) = (srho / 3) * (i ‐ 15) + mrho

nCount = WorksheetFunction.Frequency(rho, y)

For i = 2 To 29

Cells(i, 1).Value = y(i) ‐ (srho / 6)

Cells(i, 2).Value = nCount(i, 1)

End Sub

80 一つの級の幅が(標準偏差/3) = (srho/3)である。yには階級の上限値が記録されているので、そこから半幅である(srho/6)を減じることで、階級の中心値を計算している。ヒストグラムは平均値を中心に±5× 標準偏差の範囲を30階級に区切って描く。一つの級の幅が(標準偏差/3) = (srho/3)となる。

(81)

気を付けたいこと

• 出力の平均値（2.5 g∙cm

–3

）は、入力の平均値を式に代入して

得た値（2.4 g∙cm

–3

_{）と違う。}

• 95 %の信頼の水準の区間の取り方はいくらでもある。（0 ~

95 %点でも、2.5 ~ 97.5 %点でも、5 ~ 100 %点でも、信頼の水

準は95 %。）

81



_g/cm

3





出現回数

(82)

平均値の違い

• 多くの場合、モンテカルロ計算で求めた出力

の平均値と、入力の平均値をモデル式に当

てはめた値は異なる。

• JCGM 101では、モンテカルロ計算で求めた出

力の平均値を、測定量の値とすることを定め

ることを基本としている。

82

(83)

包含区間の選択

• ２章で例に見せた2.5 %点と97.5 %点を両端に

持つ区間も95 %の信頼の水準の区間である。

• JCGM 101では「最短包含区間」を提唱してい

る。これは、無数の包含区間の中で最短のも

のを選ぶということである。

• ここでは、0 %～95 %点の区間から、5 %～

100 %点の区間を、0.1 %刻みで変化させなが

ら、最短となる区間を探索する方法を紹介す

る。（他のやり方もありうる。）

83

(84)

Sub Density1Coverage()

Dim i As Long

Dim mu1 As Double

Dim sigma1 As Double

Dim x1(1 To 1000000) As Double

Dim mu2 As Double

Dim sigma2 As Double

Dim x2(1 To 1000000) As Double

Dim p As Double

Dim rho(1 To 1000000) As Double

Dim rhoi As Double

Dim rhoe As Double

Dim d As Double

Dim dmin As Double

Dim rhomin As Double

Dim rhomax As Double

（次スライドに続く）

84

(85)

mu1 = 1#

sigma1 = 0.1

For i = 1 To 1000000

x1(i) = WorksheetFunction.Norm_Inv(Rnd, mu1, sigma1)

mu2 = 10#

sigma2 = 0.1

For i = 1 To 1000000

x2(i) = WorksheetFunction.Norm_Inv(Rnd, mu2, sigma2)

p = WorksheetFunction.pi()

For i = 1 To 1000000

rho(i) = (3 / (4 * p)) * x2(i) / x1(i) ^ 3

(次スライドに続く)

85 ここまでは同じ。

(86)

rhomin = WorksheetFunction.Percentile(rho, 0#)

rhomax = WorksheetFunction.Percentile(rho, 0.95)

dmin = rhomax ‐ rhomin

For i = 1 To 50

di = 0.001 * i

rhoi = WorksheetFunction.Percentile(rho, 0# + di)

rhoe = WorksheetFunction.Percentile(rho, 0.95 + di)

d = rhoe – rhoi

If (d < dmin) Then

dmin = d

rhomin = rhoi

rhomax = rhoe

End If

Cells(1, 3).Value = rhomin

Cells(2, 3).Value = rhomax

End Sub

86 最後にrhomin、rhomaxに記録ｓれているのが最短区間の上限値と下限値である。最初の候補として、(0, 95) %点からなる区間を考える。その下限と上限を rhomin、rhomaxに入れる。dminはその区間の幅。 0.1 %ずつ区間をずらしながら、包含区間の下限値rhoiと上限値rhoe、ならびにその区間の幅dを計算している。もし、dが候補となっている最短区間の幅dminより小さかったら、新しい区間を最短区間の候補とし、rhomin、 rhomax、dminを更新する。