自由度 9 の t 分布から、乱数を 1000 個発生し、上 側の 2.5 % 点を求めるプログラム「 T_Percent 」を
3. タイプ A の不確かさの計算が違う。
71
• JCGM 101
では、 ベイズ統計学 という学
問体系に基づいて、不確かさを計算すること が基本になっている。
• JCGM 100
の従来的な不確かさ評価手法とは
タイプＡの不確かさに対する考え方が違う。
• t
分布がタイプ
Aの不確かさに現れるので、
ちょっと面倒に思える一方で、「有効自由度」
を計算する必要がなくなる。
ある測定値が正規分布に従っている。その平 均
についても、標準偏差
についても事前の 情報はないものとする。４回の測定の結果は、
9
、
12、
11、
10であった。平均
の値の事後分 布はどのように与えられるか。
=_____
= ______ ?
?
72
正規分布の平均の推定：
正規分布から得た平均の事後確率密度は、繰 り返し数ー１の自由度の
t分布の値に以下をか けて、平均値を足した値の確率密度である。
残差平方和
繰り返し数×
（繰り返し数－１）
フリークエンティスト 統計における平均値 の標準偏差
73
10.5
平均値の標準偏差を計算すると0.65となる。
JCGM 101では、ベイズ統計学に基づき、
は
「自由度４のｔ分布」の値に0.65をかけ、平 均値10.5を足した分布を持つと与える。
の事後確率密度
74
詳細は、
https://staff.aist.go.jp/k.shir ono/downloadBayes.html を参照
タイプＡの不確かさ評価
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
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)を足している。
平均と標準偏差はここではワークシー ト関数を用いて計算する。
77
ある球の密度を知るのに、直径の測定と質量の測定を する。半径の値
rは平均
1.0 cm、標準偏差
0.1 cmの正 規分布に従って測定されたとする。質量の値
mは平均
10.0 g、標準偏差
0.1 gの正規分布に従うとする。密度 の値、標準不確かさと
95 %の信頼の水準の区間を求 めよ。（密度の値は以下の式で値は
2.4 g∙cm−3。 ）
3 4 3
3 3
4 r
m r
m
g∙cm–3
分布の様子を調べるために、ヒストグラムを作成する。
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
（次スライドに続く）
78x1に半径の値を入れる。mu1、sigma1 はそれぞれ正規分布の平均と標準偏 差。
x2に質量の値を入れる。mu2、sigma2 はそれぞれ正規分布の平均と標準偏 差。
rhoに密度の値を入れる。mrhoとsrho は結果として得られる分布の平均と標 準偏差。
yはヒストグラムの級の上限値。
nCountにヒストグラムのデータを出力
する。
mu1 = 1#
sigma1 = 0.1
For i = 1 To 1000000
x1(i) = WorksheetFunction.Norm_Inv(Rnd, mu1, sigma1) Next
mu2 = 10#
sigma2 = 0.1
For i = 1 To 1000000
x2(i) = WorksheetFunction.Norm_Inv(Rnd, mu2, sigma2) Next
p = WorksheetFunction.pi() For i = 1 To 1000000
rho(i) = (3 / (4 * p)) * x2(i) / x1(i) ^ 3 Next
（次スライドに続く）
79rhoに密度の値を入れる。
x1に半径の値を入れる。mu1、sigma1 はそれぞれ正規分布の平均と標準偏 差。
x2に質量の値を入れる。
pに = 3.14…を記録する。
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
80
一つの級の幅が(標準偏差/3) = (srho/3)である。yには階級の上限値 が記録されているので、そこから半幅 である(srho/6)を減じることで、階級の 中心値を計算している。
ヒストグラムは平均値を中心に±5× 標準偏差の範囲を30階級に区切って 描く。一つの級の幅が(標準偏差/3) = (srho/3)となる。
気を付けたいこと
•
出力の平均値（
2.5 g∙cm–3）は、入力の平均値を式に代入して 得た値（
2.4 g∙cm–3）と違う。
• 95 %
の信頼の水準の区間の取り方はいくらでもある。（
0 ~ 95 %点でも、
2.5 ~ 97.5 %点でも、
5 ~ 100 %点でも、信頼の水 準は
95 %。）
g/cm3
81
出現回数
平均値の違い
•
多くの場合、モンテカルロ計算で求めた出力 の平均値と、入力の平均値をモデル式に当 てはめた値は異なる。
• JCGM 101
では、モンテカルロ計算で求めた出
力の平均値を、測定量の値とすることを定め ることを基本としている。
82
包含区間の選択
•
２章で例に見せた
2.5 %点と
97.5 %点を両端に 持つ区間も
95 %の信頼の水準の区間である。
• JCGM 101
では「最短包含区間」を提唱してい
る。これは、無数の包含区間の中で最短のも のを選ぶということである。
•
ここでは、
0 %～
95 %点の区間から、
5 %～
100 %
点の区間を、
0.1 %刻みで変化させなが ら、最短となる区間を探索する方法を紹介す る。（他のやり方もありうる。）
83
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
mu1 = 1#
sigma1 = 0.1
For i = 1 To 1000000
x1(i) = WorksheetFunction.Norm_Inv(Rnd, mu1, sigma1) Next
mu2 = 10#
sigma2 = 0.1
For i = 1 To 1000000
x2(i) = WorksheetFunction.Norm_Inv(Rnd, mu2, sigma2) Next
p = WorksheetFunction.pi() For i = 1 To 1000000
rho(i) = (3 / (4 * p)) * x2(i) / x1(i) ^ 3 Next
(
次スライドに続く
)85
ここまでは同じ。
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
Next
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を更新する。
最短区間
•
最短区間は（
1.2, 4.2）
g∙cm–3。
2.5 %点～
97.5 %点である
（
1.4, 4.6) g∙cm–3よりも、わずかに短い。
g/cm3
87
出現回数
演習
88
89
ある球の密度を知るのに、直径の測定と質量の測定 をする。半径の値
rは平均
1.0 cm、標準偏差
0.1 cmの 正規分布に従って測定されたとする。質量を
5回測定 し、
(10.3, 10.1, 10.0, 9.9, 9.7) gという結果を得た。質量 の値
mの決定に繰り返し以外の不確かさはないものと する。密度の値、標準不確かさ、
95 %の信頼の水準の 最短区間を求め、ヒストグラムを作成せよ。
3 4 3
3 3
4 r
m r
m
g∙cm–3
Sub Density_Enshu() Dim i As Long
Dim mu1 As Double Dim sigma1 As Double
Dim x1(1 To 1000000) 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 1000000) As Double Dim pi As Double
Dim rho(1 To 1000000) As Double (
次スライドに続く
)90
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
mu1 = 1#
sigma1 = 0.1
For i = 1 To 1000000
x1(i) = WorksheetFunction.Norm_Inv(Rnd, mu1, sigma1) Next
（以下省略）
End Sub
92
4.
その他・
まとめ
93
プログラミング環境①
•
エクセルＶＢＡの⻑所
•
データの管理をエクセルできる。
•
操作性という観点で親しみやすい。
•
エクセルＶＢＡの短所
•
コーディングが面倒である。
•
計算が遅い。
個人的な好みの問題が大きいと思う。他のプロ
グラミング環境も紹介しておく。
プログラミング環境②
エクセル以外のプログラミング環境
(1)：
• R：統計計算とグラフィックスのための無料の
ソフトウェア環境。
• The R Project for Statistical Computing (https://www.r‐project.org/)
•
ダウンロードサイトへのリンクあり。
•
日本語の解説も充実している。
プログラミング環境③
エクセル
VBA以外のプログラミング環境
(2)：
• NIST Uncertainty Machine：米国の国家計量機関
（NIST）が提供している無料のモンテカルロ法 に特化した計算環境。（Javaを使用）
• NIST Uncertainty Machine (https://uncertainty.nist.gov/)
•