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

Transcription

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

2 1. モンテカルロ計算で不確かさ評価 ( 原理 ) 2

3 ある球の密度を知るのに 直径の測定と質量の測定をする 半径の値 r は平均 1.0 cm 標準偏差 0.1 cm の正規分布に従って測定されたとする 質量の値 m は平均 10.0 g 標準偏差 0.1 g の正規分布に従うとする 密度の値は以下の式で値は 2.4 g cm 3 となる 標準不確かさと 95 % の信頼の水準の区間を求めよ m 3 m r r 3 g cm 3 3

4 不確かさの伝播則を用いて標準不確かさ評価する 感度係数 をrで微分すると g cm 3 4 r m g cm 4 を m で微分すると r cm -3 4

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 これらを乗じて g cm 3 と計算される の標準不確かさ u( ) u g cm 3 = 0.72 g cm 3 5

6 の拡張不確かさ U 中心極限定理から正規分布として 包含係数 2 を用いると U = 2 u( ) = 1.4 g cm 3 95 % の信頼の水準の区間 ( U, + U) = ( , ) g cm 3 = (1.0, 3.8) g cm 3 6

7 バジェットシートは以下のようになる 不確かさ要因 記号 標準不確かさ u(x i ) 半径の測定 u(r) 0.1 cm 質量の測定 u(m) 0.1 g 合成標準不確かさ拡張不確かさ 感度係数 c i 測定量の標準不確かさ c i u(x i ) (g cm 3 ) 7.2 g cm cm u c ( ) 0.72 備考 U 1.44 k = 2 7

8 モンテカルロ法でやってみる エクセルで標準偏差 0 分散 1 の正規分布はセルに =NORM.INV(RAND(),0,1) と打ち込むと得られる

9 平均 1.0 cm 標準偏差 0.1 cm の正規分布を得るには =NORM.INV(RAND(),1.0,0.1) と打ち込む

10 同じ数式を A2 から A100 までコピーします

11 平均 10.0 g 標準偏差 0.1 g の正規分布を得るには =NORM.INV(RAND(),10.0,0.1) と打ち込む

12 同じ数式を B2 から B100 までコピーします

13 A1 に入っている 半径 r B1 に入っている 質量 m から 密度 r を算出するために C1 に =3/(4*PI())*B1/A1^3 と入力します

14 同じ数式を C2 から C100 までコピーします

15 標準偏差を =STDEV.S(C1:C100) で計算します

16 95 % の信頼の水準を持つ区間の最小値は 全体の 2.5 パーセント点である PERCENTILE 関数により パーセント点は計算できる この場合 =PERCENTILE(C1:C100,0.025) とする 16

17 95 % の信頼の水準を持つ区間の最大値は 全体の 97.5 パーセント点である PERCENTILE 関数により パーセント点は計算できる この場合 =PERCENTILE(C1:C100,0.975) とする 17

18 このデータから 標準不確かさは0.83 g cm 3 95 % の信頼の水準の区間は (1.4, 4.8) g cm 3 と分かる 感度係数を使って計算した 95 % の信頼の水準の区間は (1.0, 3.8) g cm 3 18

19 さらに計算を続けると以下のヒストグラムを得る ヒストグラム : ある定められた区間にあてはまる値が何回現れたかを示す棒グラフ 観測回数 u( ) = 0.83 g cm % 点と 97.5 % 点 (1.4, 4.6) g cm 3 / g cm 3 19

20 モンテカルロ計算では 1. 感度係数は計算しない ( 感度係数が計算できないときにモンテカルロ計算は便利 ) 2. 中心極限定理は使わない ( 中心極限定理が使えないときにモンテカルロ計算は便利 ) 3. タイプ A の不確かさの計算が違う ( 後ほど )

21 便利な点 非正規性や非線形性が問題になる場合とならない場合の場合分けが必要でない 問題点 ソフトウェアの正しさ(ISO a) の証明をどのようにするのかが一般的でない モンテカルロ法による計算結果を校正結果として受け取った場合の取り扱いがどうしてよいか分からない 21

22 JCGM Evaluation of measurement data Supplement 1 to the "Guide to the expression of uncertainty in measurement" Propagation of distributions using a Monte Carlo method 22

23 2. VBA プログラミング入門 23

24 開発 タブを追加する 1 1. ファイル タブをクリック 2. オプション をクリック 3. リボンのユーザー設定 をクリック 24

25 開発 タブを追加する 2 4. メインタブ のチェックボックスの中から 開発 をチェック 5. OK ボタンをクリック 6. 開発 タブが追加された 25

26 マクロを使う 1 1. 開発 タブをクリック 2. マクロ をクリック 3. マクロ名 に適当な名称を入力し 作成 ボタンをクリックする 26

27 マクロを使う 2 4. マクロ編集画面が立ち上がる 6. 編集画面の右上の最小化ボタンをクリックして 編集画面を最小化し エクセルを表示する 5.Sub の行と End Sub の間に Msgbox( モンテカルロ法始めました ) と入力 27

28 マクロを使う 3 7. マクロ をクリック 8. 作成したマクロを選択し 実行 ボタンをクリックする 8. メッセージが表示される 9. OK ボタンをクリック 28

29 宣言 型 Sub Sengen () Dim ias Long Dim x As Double i= 1 x = 1# MsgBox (i) MsgBox (x) End Sub i を整数型として扱うことを宣言する x を実数型として扱うことを宣言する i に 1 という値を記録する 1# は整数としての 1 ではなく 実数としての のことを意味している MsgBox() はメッセージボックスに値を表示させる関数 29

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 を代入する 以下同様 30

31 データの入出力 1 Sub NyuShutsuRyoku1() Dim x As Double x = Range("A1").Value Range("B1").Value = x End Sub セルの A1 に入っているデータを x に入力する x に入っているデータをセルの B1 に出力する 31

32 データの入出力 2 Sub NyuShutsuRyoku2() Dim x As Double x = Cells(2,1).Value Cells(2,2).Value = x End Sub 上から 2 段目 左から 1 行目のセル すなわち A2 のセルからデータを x に入力する x に入っているデータを上から 2 段目 左から 2 行目 すなわち B2 のセルに出力する 32

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 で除する 33

34 簡単な関数 1 Sub KansuCos() Dim x As Double Dim y As Double x = 0 y = cos(x) MsgBox (y) End Sub cos() は cos を返す関数 簡単な関数は VBA で定義されている ( 次スライド参照 ) 34

35 簡単な関数 2 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 10 (x) は Log(x)/Log(10#) その他にもたくさんの関数が用意されている 35

36 ワークシート関数 1 エクセルで使える関数が VBAではそのままでは使えないことがある 例えば 最大値を呼び出す Max 関数は VBAでは準備されていない しかし WorksheetFunction.Max() とすることで エクセルを用いたVBAではMax 関数を使うことができる 同じように 使用できる関数が他にもいくつかある ただし ワークシート上で. が入っている関数は それを _ に置き換える必要がある場合もある ( 具体例はのちほど ) 36

37 ワークシート関数 2 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 に記録する このように配列を入力にできる関数もある 37

38 繰返し (For 文 ) Sub Kurikaeshi() Dim ias Long Dim n As Long n = 0 For i= 1 To 10 n = n + i Next 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 に更新される 38

39 演習 100 から 1000 までの整数の和はいくつになるか? これを計算し 出力する Goukei というプログラムを作れ 39

40 条件分岐 (If 文 )1 Sub JoukenBunki() Dim ias Long i= 1 If (i = 1) Then MsgBox ("iは1 ") Else MsgBox ("iは1でない ") End If End Sub If(( 条件式 )) Then ( 演算 1) Else ( 演算 2) End If の構文で条件分岐を行う ( ここで使われた条件式の i= 1 は i が 1 に等しい という意味 ) 条件式が真ならば Then と Else の間の演算を行う 条件式が偽のとき Else と End If の間の演算を行う もし 偽のときには何の演算もしないならば Else を省略して If(( 条件式 )) Then ( 演算 1) End If としてもよい 40

41 条件分岐 (If 文 )2 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 の代わりに数字を直接入力してもよい 41

42 乱数の呼び出し Sub Ransu() Dim x As Double x = Rnd MsgBox (x) End Sub Rnd は 0 から 1 の範囲の一様分布を呼び出す 42

43 一様分布 1 不確かさの計算では 最小値と最大値しか分かっていない値によく用いられる a (a + b)/2 b x

44 一様分布 2 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 の範囲での一様分布となる 44

45 ヒストグラム 1 ヒストグラムを作成するための Frequency 関数を紹介する まずはワークシートで使ってみる 1. データをある A 列 ( 下の場合 A1 から A6) に入れる 2. データ区切りの上限となる数値を B 列 ( 下の場合 B1 から B3) に小さい順に入れる 3.C 列のセルを B 列で使ったセルよりも一つ多く選択する ( 下の場合 C1 から C4) 45

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

47 ヒストグラム 4 6.C 列のデータを選択し 挿入 タブを選択 グラフ から 2-D 縦棒 の 集合縦棒 を選びます 7. グラフが表示される ( 実際にこの機能を使うときに 横軸の値は実際の測定値とは関係のない値になっているので 適切に修正する 修正の仕方はのちほど ) 47

48 ヒストグラム 5 次に VBA 中で 1000 個の 0 から 1 を範囲とする一様分布からの乱数がある場合に ヒストグラムを描く 0 1 x 48

49 ヒストグラム 6 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 Next ( 次スライドに続く ) Variant は任意の型ということ 計算の中で自動的に型が与えられる Frequency 関数の出力は Variant 型で指定しないといけない 49

50 ヒストグラム 7 ここでは ncount(1,1) が最小値から y(1) の間の配列 x 中のデータ数 ncount(2,1) が y(1) から y(2) までのデータ数 以下同様で ncount(11,1) が y(10) から最大値までのデータ数を返す For i= 1 To 10 y(i) = i* 0.1 Next ncount = WorksheetFunction.Frequency(x, y) For i= 1 To 10 Cells(i, 1) = y(i) 0.05 Cells(i, 2) = ncount(i, 1) Next End Sub 上限値ではなく 範囲中の平均値が出力されるように 0.05 を引いている ncount(1,1) から ncount(11,1) まで値が格納されているが ncount(11,1) はゼロなので 出力しないようにした 50

51 ヒストグラム 8 ヒストグラムを描いてみよう 1. 出力された B 列のデータを選択し 挿入 タブを選択 グラフ から 2 -D 縦棒 の 集合縦棒 を選ぶ この時点では横軸がデータと合っていない 2. 挿入 タブの データの選択 を選び 51

52 ヒストグラム 9 3. 立ち上がったデータソースの選択から 横軸 ( 項目 ) ラベル の 編集 ボタンをクリック 立ち上がった データ範囲の選択と書いている左の ボタンをクリック 4. この状態で A1 から A10 を選択すると 下のように選択したセルが入力される ここで 入力枠の右の をクリック 52

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

54 正規分布 1 不確かさの計算では 校正結果などの 中心極限定理がよくあてはまる値に用いられる ( 平均 ) ( 標準偏差 ) 2 は分散 54

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

56 正規分布 3 Sub SeikiBunpu() Dim mu As Double Dim sigma As Double Dim x As Double ここでmuは正規分布の平均 mu = 100# sigma = 10# x = WorksheetFunction.Norm_Inv(Rnd, mu, sigma) MsgBox (x) End Sub ここで sigma は正規分布の標準偏差 Norm.Inv 関数はワークシート上で累積分布関数の逆関数となっている ここでは WorksheetFunction とした上で. を _ に置き換える必要がある () の中に 0 ~ 1 の一様乱数 平均 標準偏差の順に入力 56

57 正規分布 4 これで計算は動く しかし ときに Rnd の値が極めて 0 あるいは 1 に近いときに WorksheetFunction.Norm_Inv(Rnd, mu, sigma) でうまく計算が行われないことがある これを避けるには Rnd が極端に 0 や 1 に近い値になるのを防ぐために Rnd の代わりに (1 2*1.E 16)*Rnd+1E 16 などとすることがある 1E 16 は小さい値で計算に支障のない値なら 他の値でもよい このテキストでは コードの簡明さを重視してこのテクニックは使用しない 57

58 演習 標準正規分布から乱数を 1000 個発生し そのヒストグラムを描くための数値を出力するプログラム Seiki_Hist を作成せよ 5 から 5 の範囲を 20 個に区切るものとする 58

59 t 分布 1 モンテカルロ計算でタイプ A の不確かさ評価をするときに使う分布 ( 通常の不確かさ評価とは使い方が違うので注意 詳しくは後述 ) = 1 = 3 = 5 は自由度と呼ばれる変数 これを決めると t 分布の形が決まる 59

60 t 分布 3 Sub TBunpu() Dim nu As Long Dim x As Double nuは自由度 nu = 4 x = WorksheetFunction.T_Inv(Rnd, nu) MsgBox (x) End Sub ここでも 逆関数法を使う T.Inv 関数はワークシート上で累積分布関数の逆関数となっている () の中に 0 ~ 1 の一様乱数 自由度の順に入力 60

61 三角分布 1 不確かさ評価では 三角分布の情報が事前に与えられているときに使う 多くの場合には 同一の一様分布を 2 つ重ね合わせた不確かさ要因である 61

62 三角分布 2 独立に 0 ~ 1 の範囲の一様分布に従う 2 つの値を x y とすると z = x y は 1 ~ 1 の範囲の三角分布である 0 x 1 1 z = x y 1 0 y 1 62

63 三角分布 3 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 を足す 63

64 U 字分布 1 不確かさ評価では sin 関数に従って規則的に時間変動する変数に対して用いる b x a 0 時間 a x b 64

65 U 字分布 2 Sub UjiBunpu() Dim a As Double Dim b As Double Dim x As Double a = 100 b = 100 p = WorksheetFunciton.Pi() x = (b a)/2*sin(2*p*rnd)+(b+a)/2 MsgBox (x) End Sub ここで a と b は U 字分布の最小値と最大値 p には = を代入しておく p に を代入した上で (b a)/2*sin(2*p*rnd) とすると 平均 0 全幅 (b a) の U 字分布となる それに平均値の (b+a)/2 を足す 65

66 パーセンタイル点 1 標準正規分布からの 1000 個のサンプルを導き 上側 2.5 パーセント点を計算する? = 1 = % 66

67 パーセンタイル点 2 Sub SeikiBunpuPercentile() Dim ias 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) Next y = WorksheetFunction.Percentile(x, 0.975) MsgBox (y) End Sub Percentile 関数に ( 配列 確率 ) を入力することで 下側のパーセント点を求めることができる 上側 100a % 点は下 側 ( a) % パーセントである 67

68 演習 自由度 9 の t 分布から 乱数を 1000 個発生し 上側の 2.5 % 点を求めるプログラム T_Percent を作成せよ 68

69 3. モンテカルロ計算で不確かさ評価 ( 応用 ) 69

70 モンテカルロ計算では 1. 感度係数は計算しない ( 感度係数が計算できないときにモンテカルロ計算は便利 ) 2. 中心極限定理は使わない ( 中心極限定理が使えないときにモンテカルロ計算は便利 ) 3. タイプ A の不確かさの計算が違う

71 JCGM 101では ベイズ統計学という学問体系に基づいて 不確かさを計算することが基本になっている JCGM 100の従来的な不確かさ評価手法とはタイプAの不確かさに対する考え方が違う t 分布がタイプAの不確かさに現れるので ちょっと面倒に思える一方で 有効自由度 を計算する必要がなくなる 71

72 ある測定値が正規分布に従っている その平均 についても 標準偏差 についても事前の情報はないものとする 4 回の測定の結果は であった 平均 の値の事後分布はどのように与えられるか =? =? 72

73 正規分布の平均の推定 : 正規分布から得た平均の事後確率密度は 繰り返し数ー 1 の自由度の t 分布の値に以下をかけて 平均値を足した値の確率密度である 残差平方和 フリークエンティスト統計における平均値の標準偏差 繰り返し数 ( 繰り返し数 -1) 73

74 平均値の標準偏差を計算すると 0.65 となる JCGM 101 では ベイズ統計学に基づき は 自由度 4 の t 分布 の値に 0.65 をかけ 平均値 10.5 を足した分布を持つと与える の事後確率密度 詳細は ono/downloadbayes.html を参照

75 タイプ A の不確かさ評価 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 にモンテカルロ法で得たサンプルを記録する 75

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) Next MsgBox (WorksheetFunction.StDev_S(x1) ) End Sub 自由度 4(= n1 1) の t 分布の値に計算した平均値の標準偏差の値 (s1) をかけて 平均値 (m) を足している 76

77 ある球の密度を知るのに 直径の測定と質量の測定をする 半径の値 r は平均 1.0 cm 標準偏差 0.1 cm の正規分布に従って測定されたとする 質量の値 m は平均 10.0 g 標準偏差 0.1 g の正規分布に従うとする 密度の値 標準不確かさと 95 % の信頼の水準の区間を求めよ ( 密度の値は以下の式で値は 2.4 g cm 3 ) m 3 m r r 3 g cm 3 77

78 分布の様子を調べるために ヒストグラムを作成する Sub Density1Histgram() Dim ias Long Dim mu1 As Double Dim sigma1 As Double Dim x1(1 To ) As Double Dim mu2 As Double Dim sigma2 As Double Dim x2(1 To ) As Double Dim p As Double Dim rho(1 To ) As Double Dim mrho As Double Dim srho As Double Dim ncount As Variant Dim y(1 To 29) As Double する ( 次スライドに続く ) 78 x1 に半径の値を入れる mu1 sigma1 はそれぞれ正規分布の平均と標準偏差 x2 に質量の値を入れる mu2 sigma2 はそれぞれ正規分布の平均と標準偏差 rho に密度の値を入れる mrho と srho は結果として得られる分布の平均と標準偏差 y はヒストグラムの級の上限値 ncount にヒストグラムのデータを出力

79 mu1 = 1# sigma1 = 0.1 For i= 1 To x1(i) = WorksheetFunction.Norm_Inv(Rnd, mu1, sigma1) Next mu2 = 10# sigma2 = 0.1 For i= 1 To x2に質量の値を入れる x2(i) = WorksheetFunction.Norm_Inv(Rnd, mu2, sigma2) Next pに = 3.14 を記録する p = WorksheetFunction.pi() For i= 1 To rhoに密度の値を入れる rho(i) = (3 / (4 * p)) * x2(i) / x1(i) ^ 3 Next ( 次スライドに続く ) 79 x1 に半径の値を入れる mu1 sigma1 はそれぞれ正規分布の平均と標準偏差

80 mrho = WorksheetFunction.Average(rho) srho = WorksheetFunction.StDev_S(rho) For i= 1 To 29 y(i) = (srho / 3) * (i 15) + mrho Next 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) Next End Sub ヒストグラムは平均値を中心に ±5 標準偏差の範囲を 30 階級に区切って描く 一つの級の幅が ( 標準偏差 /3) = (srho/3) となる 一つの級の幅が ( 標準偏差 /3) = (srho/3) である y には階級の上限値が記録されているので そこから半幅である (srho/6) を減じることで 階級の中心値を計算している 80

81 気を付けたいこと 出力の平均値 (2.5 g cm 3 ) は 入力の平均値を式に代入して得た値 (2.4 g cm 3 ) と違う 95 % の信頼の水準の区間の取り方はいくらでもある (0 ~ 95 % 点でも 2.5 ~ 97.5 % 点でも 5 ~ 100 % 点でも 信頼の水準は95 % ) 出現回数 3 g/cm 81

82 平均値の違い 多くの場合 モンテカルロ計算で求めた出力の平均値と 入力の平均値をモデル式に当てはめた値は異なる JCGM 101では モンテカルロ計算で求めた出力の平均値を 測定量の値とすることを定めることを基本としている 82

83 包含区間の選択 2 章で例に見せた 2.5 % 点と 97.5 % 点を両端に持つ区間も 95 % の信頼の水準の区間である JCGM 101 では 最短包含区間 を提唱している これは 無数の包含区間の中で最短のものを選ぶということである ここでは 0 %~95 % 点の区間から 5 %~ 100 % 点の区間を 0.1 % 刻みで変化させながら 最短となる区間を探索する方法を紹介する ( 他のやり方もありうる ) 83

84 Sub Density1Coverage() Dim ias Long Dim mu1 As Double Dim sigma1 As Double Dim x1(1 To ) As Double Dim mu2 As Double Dim sigma2 As Double Dim x2(1 To ) As Double Dim p As Double Dim rho(1 To ) 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 x1(i) = WorksheetFunction.Norm_Inv(Rnd, mu1, sigma1) Next mu2 = 10# sigma2 = 0.1 For i= 1 To x2(i) = WorksheetFunction.Norm_Inv(Rnd, mu2, sigma2) Next p = WorksheetFunction.pi() For i= 1 To rho(i) = (3 / (4 * p)) * x2(i) / x1(i) ^ 3 Next ここまでは同じ ( 次スライドに続く ) 85

86 rhomin = WorksheetFunction.Percentile(rho, 0#) rhomax = WorksheetFunction.Percentile(rho, 0.95) dmin = rhomax rhomin For i= 1 To 50 di = * i の区間の幅 rhoi = WorksheetFunction.Percentile(rho, 0# + di) rhoe = WorksheetFunction.Percentile(rho, di) d = rhoe rhoi If (d < dmin) Then dmin = d rhomin = rhoi rhomax = rhoe End If Next Cells(1, 3).Value = rhomin Cells(2, 3).Value = rhomax End Sub 最初の候補として (0, 95) % 点からなる区間を考える その下限と上限を rhomin rhomax に入れる dmin はそ 0.1 % ずつ区間をずらしながら 包含区間の下限値 rhoi と上限値 rhoe ならびにその区間の幅 d を計算している もし d が候補となっている最短区間の幅 dmin より小さかったら 新しい区間を最短区間の候補とし rhomin rhomax dmin を更新する 最後に rhomin rhomax に記録 s れているのが最短区間の上限値と下限値である 86

87 最短区間 最短区間は (1.2, 4.2) g cm % 点 ~97.5 % 点である (1.4, 4.6) g cm 3 よりも わずかに短い 出現回数 3 g/cm 87

88 演習 88

89 ある球の密度を知るのに 直径の測定と質量の測定をする 半径の値 r は平均 1.0 cm 標準偏差 0.1 cm の正規分布に従って測定されたとする 質量を 5 回測定し (10.3, 10.1, 10.0, 9.9, 9.7) g という結果を得た 質量の値 m の決定に繰り返し以外の不確かさはないものとする 密度の値 標準不確かさ 95 % の信頼の水準の最短区間を求め ヒストグラムを作成せよ m 3 m r r 3 g cm 3 89

90 Sub Density_Enshu() Dim ias Long Dim mu1 As Double Dim sigma1 As Double Dim x1(1 To ) As Double Dim n2 As Long Dim rx2(1 To 5) As Double Dim rx2m As Double Dim rx2s As Double Dim x2(1 To ) As Double Dim pi As Double Dim rho(1 To ) As Double ( 次スライドに続く ) 90

91 Dim mrho As Double Dim srho As Double Dim ncount As Variant Dim y(1 To 29) 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 ( 次スライドに続く ) 91

92 mu1 = 1# sigma1 = 0.1 For i= 1 To x1(i) = WorksheetFunction.Norm_Inv(Rnd, mu1, sigma1) Next ( 以下省略 ) End Sub 92

93 4. その他 まとめ 93

94 プログラミング環境 1 エクセル VBA の 所 データの管理をエクセルできる 操作性という観点で親しみやすい エクセル VBA の短所 コーディングが面倒である 計算が遅い 個人的な好みの問題が大きいと思う 他のプログラミング環境も紹介しておく

95 プログラミング環境 2 エクセル以外のプログラミング環境 (1): R: 統計計算とグラフィックスのための無料のソフトウェア環境 The R Project for Statistical Computing ( project.org/) ダウンロードサイトへのリンクあり 日本語の解説も充実している

96 プログラミング環境 3 エクセルVBA 以外のプログラミング環境 (2): NIST Uncertainty Machine: 米国の国家計量機関 (NIST) が提供している無料のモンテカルロ法に特化した計算環境 (Javaを使用) NIST Uncertainty Machine ( プログラミングの知識が全くない人にはこちらの方がとっつきやすいかも知れない

97 まとめ コメント モンテカルロ法による不確かさの算出方法について解説した 非線形性や非正規性が問題になる場合には 非常に有用なツールである 今回は触れられなかったが JCGM 101では通常のGUMの方式と モンテカルロ計算を比較することで GUMの方式の正当性を示すような妥当性評価方法も提案されている JCGM 101の内容などもよく確認し 不確かさ評価にモンテカルロ計算をうまく活用して欲しい