不確かさ評価の基礎
名古屋大学 山本章夫

目次
1.はじめに ................................................................................................................................ 3
2.科学における予測の不確かさ ................................................................................................. 3
3. 誤差と不確かさ .................................................................................................................... 4
3.1 誤差と不確かさの違い ......................................................................................................... 4
3.2 誤差の種類 ........................................................................................................................... 6
3.3 不確かさ評価の種類 ............................................................................................................. 8
4.確率密度関数 ....................................................................................................................... 10
4.1 確率関数 ............................................................................................................................. 11
4.2 確率密度関数 ...................................................................................................................... 11
4.3 累積分布関数 ...................................................................................................................... 13
5.期待値、分散、標準偏差 ...................................................................................................... 14
5.1 期待値................................................................................................................................. 14
5.2 分散 .................................................................................................................................... 15
5.3 標準偏差 ............................................................................................................................. 17
5.4 モーメント ......................................................................................................................... 17
6.母集団の推定 ....................................................................................................................... 17
6.1 母集団とは ......................................................................................................................... 17
6.2 平均値の不確かさ(標準偏差).............................................................................................. 18
6.3 標準偏差の推定 .................................................................................................................. 19
6.4 信頼係数と信頼区間 ........................................................................................................... 21
6.5 推定臨界下限値増倍率 ....................................................................................................... 25
6.6 注意事項 ............................................................................................................................. 25
7.共分散、相関係数、分散共分散行列..................................................................................... 27
8.不確かさの伝搬 .................................................................................................................... 291

8.1 概要 .................................................................................................................................... 29
8.2 多変数関数の誤差伝搬 ....................................................................................................... 30
8.3 二変数の場合の誤差伝搬.................................................................................................... 31
8.4 不確かさ伝搬の SANDWICH FORMULA .................................................................................. 33
9.まとめ .................................................................................................................................. 35
謝辞 ........................................................................................................................................ 35
APPENDIX. 推定臨界下限増倍率の算出法 ............................................................................. 36
A.1 必要となる統計学知識 ...................................................................................................... 36
A.2 実効増倍率の母集団の平均値と標準偏差がわかっている場合 ......................................... 37
A.3 実効増倍率の母集団の平均値と標準偏差が不明であるが、これらが十分に精度良く推定できている場合 ........................................................................................................................... 38
A.4 実効増倍率の母集団の平均値が不明であり、標準偏差のみがわかっている場合 ............ 40
A.5 実効増倍率の母集団の平均値も標準偏差も不明な場合 .................................................... 442

1.はじめに
次の質問に答えられるだろうか？
①不確かさと誤差の違いについて説明せよ②複数回測定を行って得られた測定値の平均値の不確かさについて説明せよ③不確かさには、性質の異なる二つの種類がある。これらについて説明せよ①～③全てに正確に答えられる人はあまりいないのではないだろうか。原子力は、主として設計(計算)で製品の性能を保証している。そのため、計算結果がどれぐらい信頼できるか(疑わしいか)についての定量的な評価がきわめて重要になることは明らかである。不確かさ評価は、計算結果の信頼性を定量的に保証するために不可欠な考え方である。
本稿では、不確かさ評価を行うにあたって必要となる基礎知識について解説することを試みる。

2.科学における予測の不確かさ
科学の世界では、実世界における現象は、いくつかの過程をへて「予測」される。まず、測定結果に基づき、理論モデルが作成される。これらの理論モデルを(多くの場合は近似的な)解析モデルに置き直し、この解析モデルを用いて数値解析などにより予測結果を得る。
このように、予測は複雑な手順を踏んで行われるものであり、その結果が「疑わしくなりえる」要素には以下のように様々なものがある。
①観察(測定)結果：現実世界の測定は有限回しか行うことができない。また、測定は、長さであれば定規、温度であれば温度計といった何らかの補助手段を用いて行われ数値化されるが、これらの補助手段により得られる値は、当然ながら厳密なものではない。また、観察したい物理量が直接測定できず、何らかの補正や換算係数を用いて測定値を得る場
合もある。これらのことは、出発点となる測定結果にはある程度の「疑わしさ」が存在することを示している。
②理論モデル：測定結果に基づき理論モデルを構築する過程においても、我々の知識が完全であることはあり得ず、このために理論モデルには「疑わしさ」が必ず存在する。原子に関する理論モデルが変遷してきたことは、この一例である。
③解析モデル：理論モデルは、そのままでは解くことができない場合が多い。例えば、
Boltzmann 輸送方程式はその一例である。そのため、理論モデルに基づき、解析モデルを作り、結果を得ることとなる。この過程では、様々な離散化、仮定、単純化などが行われるが、これらはいずれも予測の「疑わしさ」に結びつくものである。
④入力データ：解析を行うためには、解析モデルに与える入力データが必要になる。例え3

ば原子炉の炉心解析では、燃料の組成、形状、温度、原子炉出力、と言ったパラメータが入力データとして必要となる。これらのパラメータは、全て測定により得られたものであることから、不確かさが必ず存在し、この不確かさは予測の「疑わしさ」に結びつく。
以上のことから、程度は様々であるものの、予測値は「疑わしい」ものであることがわかる。従って、予測値を有効に活用するためには、「ではどの程度疑わしいのか」に対して定量的に答える必要がある。これが不確かさ評価である。

3. 誤差と不確かさ
3.1 誤差と不確かさの違い
誤差と不確かさは、よく混同されるが、異なるものである。
誤差(error)：測定(計算)しようとする対象についての、測定(計算)された値と真の値との差(測定(計算)値－真の値)を言う。一般に、真の値は測定によっても得ることはできない。
これは、測定にも不確かさが必ず含まれているからである。従って、定義通りの「誤差」
は、永久に求めることはできない。しかし、十分に注意を払って得た測定値は、真の値に十分近いと考えて良いであろうから、炉物理の分野では、「計算値－測定値」を一般的に計算の誤差と呼んでいる。
不確かさ(uncertainty)：簡単に言うと、測定(計算)値のばらつきを特徴付けるパラメータである。別の言い方をすると、測定(計算)値の疑わしさを数値で表したものである。より正確には、測定(計算)の結果から得られる真の値の「候補」のばらつきを特徴付けるパラメー
タである。炉物理の分野では、不確かさは、誤差の標準偏差として評価することが一般的である。
「誤あるパラメータに対する誤差 X i に対し、誤差の平均( X )、標準偏差(不確かさ) σ ( X ) 、差の平均」の標準偏差 σ ( X ) の一般的な計算式を以下にまとめる。ここで、n は計算もしく
は測定回数を示す。なお、標準偏差の計算方法の詳細については、5 章で説明を行う。
n

X =

∑X
i =1

i

(3.1)

n
n

σ(X ) =∑( X − X )
i =1

i

2

(3.2)

n −1

4

σ(X ) =σ (Xi )

(3.3)

n

■例１：ここに、ある決定論的な核特性解析コードがあるとする。5 つの異なる臨界実験について解析を行った結果、実効増倍率は 1.001, 1.002, 0.998, 1.003, 1.000 となった。それぞれの臨界実験において測定された実効増倍率が 1.000 であるとき、この核特性解析コードの誤差と不確かさはどのように評価できるか？なお、実効増倍率の測定値の信頼性は十分に高く、不確かさは十分に小さいものとする。
測定値の信頼性が高いため、測定値を真の値と見なすと、誤差はそれぞれの炉心に対して、0.1%Δk/k, 0.2%Δk/k, -0.2%Δk/k, 0.3%Δk/k, 0.0%Δk/k である。従って実効増倍率の誤差の平均値は(3.1)式より0.001 + 0.002 + ( −0.002) + 0.003 + 0.000= 0.0008
5

また標準偏差(不確かさ)は(3.2)式より(1.001 − 1.0008) 2 +(1.002 − 1.0008) 2 +(0.998 − 1.0008) 2 +(1.003 − 1.0008) 2 +(1.000 − 1.0008) 2(5 − 1)= 0.0019

この例における不確かさは、核データ、炉心解析コードのモデルや数値解析アルゴリズ
ム、実験体系や使用されている材質の不確かさｓなどに起因したものである。□■例 2：モンテカルロコードで、ある臨界実験(1 体系)の解析を行った。初期乱数を変えて 5回解析を行った結果、実効増倍率は 1.001, 1.002, 0.998, 1.003, 1.000 となった。この臨界実験における実効増倍率の測定値は 1.000 であり、十分に信頼性は高いものとする。このモンテカルロコードの誤差と不確かさはそれぞれいくらか。
誤差はそれぞれの炉心に対して、0.1%Δk/k, 0.2%Δk/k, -0.2%Δk/k, 0.3%Δk/k, 0.0%Δk/kである。従って、五回の計算における誤差の平均値は以下のように計算できる。
0.001 + 0.002 + ( −0.002) + 0.003 + 0.000= 0.00085

また、標準偏差は、(1.001 − 1.0008) 2 +(1.002 − 1.0008) 2 +(0.998 − 1.0008) 2 +(1.003 − 1.0008) 2 +(1.000 − 1.0008) 2(5 − 1)= 0.0019となる。

5この例では、一つの体系について、初期乱数を変えて計算を行っている。これら 5 回の計算から得られた誤差の標準偏差をモンテカルロコードの不確かさとすることは妥当であろうか？この例では、同じ体系を計算しているため、決定論的手法を用いれば、全く同じ
実効増倍率が得られるはずである。しかし、モンテカルロ法では、統計誤差のため、毎回異なる実効増倍率が得られる。従って、5 回の計算結果から得られた標準偏差(0.0019)は、統計誤差に基づく不確かさを与えることとなる。
せっかく初期乱数を変えて 5 回も計算したのであるから、それを有効に活用したいと思うのは自然であろう。そのためには、5 回の計算の誤差の平均をモンテカルロ法の誤差とするのが良いと考えられる。そこで次に、誤差の平均値がどの程度確からしく見積もられているか、つまり、「誤差の平均値」の不確かさがどの程度あるかについて考えてみよう。今回は、5 回の計算により平均を求めているので、（誤差の標準偏差でなく)「誤差の平均値」の標準偏差は(3.3)式より0.00195= 0.00086 である。この結果から、この炉心に関して、この モ ン テ カ ル ロ コ ー ド の 計 算 誤 差 が 1.0008 ± 0.00086 の 範 囲 に 入 る 確 率 が 約 68% 、1.0008 ± 2 × 0.00086 の範囲に入る確率が約 95%、 1.0008 ± 3 × 0.00086 の範囲に入る確率が約 99.7%となる、と書きたいところであるが、これは正しくない。正しい評価については、
6.3 節で説明する。□
ここで評価した誤差の平均値は、この炉心に対し、初期乱数を変えてモンテカルロコー
ドの計算を多数繰り返し実行することにより低減することができる。つまり、「誤差の平均値」に不確かさが含まれているのは、今回の(5 回の)計算では、(初期乱数を変えて無限回モ
ンテカルロコードの計算を繰り返した)母集団の一部をサンプリングしていることにより発生したものであると解釈できる。
3.2 誤差の種類
前節では、明示しなかったが、誤差には、偶然誤差と系統誤差の二種類がある。
偶然誤差：さいころを振ったとき出る目のように、本来、偶然に支配される要因によって発生する誤差のこと。計算に含まれる偶然誤差は、測定や計算を繰り返すことにより見積もることが可能である。なお、偶然誤差自体は、数多くの測定を行ったとしても減らすことはできない。ただし、誤差の平均値については、数多くの測定により、より正確な、言い換えるとより不確かさが小さい結果を得ることができる。放射性核種の崩壊、
原子核と中性子の相互作用などは、偶然誤差を含んでいる。モンテカルロコードにおいて、中性子が原子核と反応するときの不確かさそのものは計算を繰り返し行っても低減することはできない。一方、実効増倍率や反応率などの統計誤差も偶然誤差に起因するものであるが、統計誤差は「平均値に対する不確かさ」であることから、計算回数(つまりヒストリー数)を増やすことによりこれを低減することができる。

6

系統的誤差：毎回の計算に同じ影響を与える誤差。この誤差の要因が明らかになっているかどうかは分からない。炉心解析においては、例えば拡散近似、メッシュ誤差、エネル
ギー群の縮約誤差などが系統誤差になり得る。計算に含まれる系統的誤差を評価するためには、(十分に信頼できる)測定値との比較、他の計算方法による計算結果との比較などを行う必要がある。
また、これらの誤差の定義に関連して、正確さに関する言葉の定義を以下に示す。
正確度(accuracy)：測定(計算)値が、真値にどの程度近いのかを示す尺度。系統的誤差の小さいことを示す。
精度(precision)：測定(計算)値を複数回行ったときの値のばらつきを示す尺度。偶然誤差が小さく再現性が高いことを示す。なお、真値とのずれは考慮していない。
図 3.1 にダーツを投げたときの結果を示す。(a)～(d)の分布は以下のように解釈できる。
(a)偶然誤差小、系統誤差小。もしくは正確度高、精度高。
(b)偶然誤差小、系統誤差大。もしくは正確度低、精度高。
(c)偶然誤差大、系統誤差小。もしくは正確度高、精度低。
(d)偶然誤差大、系統誤差大。もしくは正確度低、精度低。
一方、真の値が分からない場合、図 3.2 のようになり、偶然誤差が大きいか小さいか、精度が高いか低いかは判断することができるが、系統誤差の大小、正確度の高低は判断でき
ない。
正確度：高

正確度：低

正確度：?

正確度：?

(a)

(b)

(c)

(d)

精度：高

精度：高
(a)

(b)

精度：低

精度：低
(c)

(d)

図 3.1 偶然(ランダム)誤差と系統的誤差

図 3.2 偶然(ランダム)誤差と系統的誤差(真の値が分からない場合)

7

3.3 不確かさ評価の種類
不確かさ評価の方法は、誤差の種類によって異なる。以下に不確かさを評価する二つの方法を示す。
①タイプ A 評価法：計算を複数回実施し、統計的処理によって不確かさを推定する。
②タイプ B 評価法：タイプ A 以外の手法を用いて不確かさを推定する。
以下にいくつか例を示す。
a)様々な炉心に対して決定論的手法による解析を実施し、各炉心に対する計算結果の誤差を統計処理して不確かさを評価する。(タイプ A)b)モンテカルロコードの計算を(初期乱数を変化させて)複数回実行し、その結果を統計処理して不確かさを評価する。(タイプ A)c)様々な炉心体系で拡散計算と輸送計算の結果を比較し、統計処理によって拡散計算の不確
かさを推定する。(タイプ A)d)断面積の共分散データ(不確かさ)を用いて断面積をランダムにサンプリングし、複数の断面積を用いて炉心解析を実施し、その結果を統計処理して炉心特性の不確かさを評価する。(タイプ A)e)断面積の感度係数を用いて断面積データに含まれる共分散データを伝播させ、炉心特性の不確かさを評価する。(タイプ B)
タイプ A 評価法は、統計的手法を用いることから、偶然誤差に基づく不確かさを評価するための手法と見られることがある。しかし、a)や c)に示すように、タイプ A は必ずしも偶然誤差に対する不確かさ評価を行う手法ではなく、系統誤差に対する不確かさ評価に用いられることもある。
■例 3：巻き尺を使用して測定した 1 本のひもの長さの不確かさを評価する 1。
ステップ 1：どのような不確かさを考慮しなければならないかを検討する。
①巻き尺自体の不確かさ
・補正の必要性、校正がなされているか、校正で示された不確かさの定義
・巻き尺が伸びたり縮んだりしないか(校正してからの変化はどれぐらいか)
・温度や湿度などの環境が巻き尺の長さに影響を与えないか・分解能はいくらか(ここでは 1mm と仮定)
②測定対象の不確かさ
・ひもはまっすぐ置かれているか「不確かさの入門ガイド」、ASG104，独立行政法人製品評価技術基盤機構、平成 26 年 5月, p.21．
18

・温度や湿度などの環境がひもの長さに影響を与えないか
・ひもの端は明確に分かるか③測定方法の不確かさ
・ひもの端と巻き尺の端をきちんと合わせられるか・巻き尺とひもを並行に並べられるか・同じ測定を同じ条件で繰り返すことができるかステップ 2：測定の実施。測定を 10 回行って、平均値として 5.017m、長さの標準偏差として、0.0021m が得られた。一般的には、これをひもの長さの不確かさとする場合も多いかも知れないが、これだけでは不確かさの一部でしかないことに注意。
ステップ 3：ステップ 1 で考慮した各条件の不確かさを評価する。
・巻き尺は校正済みで補正は不要。校正の不確かさ(測定誤差の標準偏差)はカタログ値から長さの 0.05%である。従って、校正の不確かさに基づく測定の不確かさは 5.017m×0.0005=0.0025m である。
・巻き尺の最小目盛りは 1mm であることから、目盛り読み取りの誤差は 0.5mm を超えることはないと推測できる。すなわち、目盛り読み取りの誤差は±0.5mm の範囲に入ると推測できる。目盛り読み取りの誤差を一様分布と見なすと、これに基づく不確かさ(標準偏差)は、一様分布の上限と下限の半値幅である 0.0005m を 3 で割ったもの、つまり 0.0005 / 3 ≅ 0.0003 m である 2。
・巻き尺はまっすぐ置かれていると推定して、不確かさは 0m とする。
・巻き尺の温度や湿度に対する変化は無視できるとして、不確かさは 0m とする。
・ひもをまっすぐに置くことは難しいため、実際の長さより短めに測定値が出やすい。
この不確かさを長さの 0.1%、つまり 0.005m と想定する。この誤差がどのように分布しているかは分からないため、これを一様分布と仮定すると3、標準偏差は0.005 / 3 ≅ 0.0029 となる。
・10 回の測定を行っているため、測定値の平均の不確かさは、 0.0021 / 10 ≅ 0.0007 m
ステップ 4：ステップ 3 で求めた不確かさが独立しているかどうかを確認する。
ここでは、これらの不確かさが独立しているとする。なお、不確かさに相関がある場合については、7 章で取り扱う。
ステップ 5：既知の補正を全て含んだ測定値を評価する。
一様分布の標準偏差の計算については、例 15 を参照のこと。
短い側に測定値が出る可能性が高いと思われるが、不確かさとして、±0.005m を見込み、上限値と下限値の半値幅としては 0.005m を仮定する。

239

今回は、ひもの曲がりに起因する補正(系統誤差)0.005m を考慮して、5.017+0.005m=5.022mとなる。
ステップ 6：全ての不確かさを合成する。
今回は、それぞれの不確かさに相関がないとしているので、合成された不確かさ 4は0.00252 + 0.00032 + 0.0029 2 + 0.0007 2 = 0.0039
ステップ 7：95%の信頼水準のもとで、ひもの長さの不確かさは 0.0039×2=0.0078m となる
5。

ステップ 8：以上の結果をまとめると、ひもの長さについて、次のような記載となる。
ひもの長さは 5.022 ± 0.0078 m である。長さの不確かさは標準不確かさを 2 倍したものであり、信頼水準は 95%となる。この長さは、水平に置かれたひもを 10 回測定した値の平均値である。この結果は、測定時にひもが完全にはまっすぐに置かれていなかった影響を補正している。不確かさは「測定の不確かさに関する入門ガイド」による方法により測定した。
表 3.1 不確かさ計算のまとめ
不確かさの要因

値(±)

確率分布

除数

標準不確かさ

校正の不確かさ

0.005

正規

2

0.0025

分解能(巻き尺の最小目盛)

0.0005一様√30.0003

ひもの置き方

0.005

一様

√3

0.0029

10回繰り返し測定の平均値の標準不確かさ

0.0007正規

10.0007

合成標準不確かさ

正規と仮定

0.003995%信頼水準での不確かさ

正規と仮定

0.0078

なお、上記の結果は、表 3.1 のような形でまとめることができる。□

4.確率密度関数
不確かさを評価する際、対象とするパラメータがどのような確率で分布しているかはき
わめて重要な情報である。以下では、対象とするパラメータが飛び飛びの値を取る(離散的な)場合と連続量である場合に分けて、確率の分布を表す関数について説明する。

4
5

不確かさの伝搬については、第 8 章参照。
95%については、例 8 を参照のこと。また、信頼水準については、6.3 節を参照。
10

4.1 確率関数
対象とするパラメータが離散的な場合、それぞれのパラメータを取る確率を表す関数を確率関数 6と呼ぶ。確率関数 p ( xi ) は以下の性質を持つ。
①あるパラメータを取る確率は 0～1 になるため、 1 ≥ p ( xi ) ≥ 0②必ずどれかのパラメータを取るとることから、∑ p( x ) = 1
i

i

(4.1)
(4.2)

■例 4：一つのさいころを振った。このさいころがいわゆる「グラサイ」7でなければ、1～
6 の目が出る確率は等しく 1/6 である。従って、サイコロの目に対する確率関数は、

1p( xi ) = , xi = {1,2,3,4,5,6}6として与えられる。この確率関数は①、②の性質を共に満たしている。□■例 5：さいころを振って偶数が出れば+1，奇数が出れば-1 だけ x 軸上を移動することを考える。原点から出発し、1 回試行したときに位置 xi をとる確率関数は1, xi = {−1,1}2p( xi ) = 0, xi = {0}p ( xi ) =2 回試行したときには、

1, xi = {−2,2}4p( xi ) = 0, xi = {−1,1}p ( xi ) =p ( xi ) =1, xi = {0}2N(偶数)回試行したときには、,いずれの試行回数の場合についても、確率関数は(4.1)、(4.2)の性質を共に満たしている。□
4.2 確率密度関数
対象とするパラメータが連続量の場合、あるパラメータを取る確率を表す関数を確率密度関数とよぶ。確率密度関数 p (x ) は、以下の性質を持つ。

67

離散的な変数に対する「確率密度関数」は確率関数、または確率質量関数と呼ばれる。
特定の目が出やすくなっているいかさま用のさいころ。
11

①単位入力値あたりの確率として与えられる。これは、重さの密度が単位体積あたりの重量として定義されているのと同じである。
②正の値を取る。つまり、 p ( x ) ≥ 0(4.3)③入力値の全範囲で積分すると、1 になる。つまり、④ a < x < b を取る確率は、 P ( a < x < b) =∞∫ p( x )dx = 1−∞b∫ p( x )dx
a

(4.4)
(4.5)

■例 6：10 分ごとに電車が来る駅に行ったときの待ち時間の確率密度関数は、待ち時間を分で表したときに1, 0 ≤ x < 1010p( x ) = 0, 10 ≤ x
p( x ) =となる。この確率密度関数は、上記の(4.3), (4.4), (4.5)の性質を満たしている。□
さて、連続量の変数に対する確率密度関数を理解する上でやっかいなのが、確率密度関数の値が確率そのものを表しているのではないことである。これは④の性質に関連している。この点について理解するため、次の問いを考えてみよう。
「密度 10g/cm3 の物質がある。
この物質の重さはいくらか？」
この問いに答えるためには、物質の体積が必ず必要である。確率密度関数→重さの密度、
確率→重さのように関連づけて考えると、体積に対応する量を与えないと確率が計算できないことは明らかである。上記の問題で確率密度関数の値が 1/10 であることは、例えば待ち時間が 3 分 30 秒である確率が 1/10 であることを示しているのではない。1/10 に待ち時間の範囲を乗じることにより、初めて確率を計算することができる。例えば、待ち時間が 3
分 30 秒から 4 分の間の 1/2 分間である確率は、

1 1 1
× =
= 5%
10 2 20であるし、待ち時間が 5 分から 10 分の間の 5 分間である確率は

1
1× 5 = = 50%
10
2
である。
■例 7：巨視的断面積Σt の媒質に x=0 から中性子が入射した。 x ~ x + dx でこの中性子が媒質と初めて衝突する確率は、以下のように求められる。
まず、位置 x まで衝突しない確率は、exp( −Σ t x ) である。さらに、位置 x まで衝突せずに飛んできた中性子が x ~ x + dx で媒質と衝突する確率は、dx がきわめて小さいとして Σ t dx

12

で 与 え ら れ る 。 こ れ か ら 、 x ~ x + dx で 中 性 子 が 初 め て 媒 質 と 衝 突 す る 確 率 は 、exp( −Σ t x )Σ t dx である。従って、中性子の衝突に関する(単位飛行距離あたりの)確率密度関数は、

p( x ) = Σ t exp( −Σ t x )(4.6)

である。また、x を 0～∞まで積分すると、∞∫ Σ exp( −Σ x )dx = 1
0

t

(4.7)

t

である。この確率密度関数は、①、②、③の性質を満たしている。なお、例 6 と同じように、位置 x で中性子が初めて媒質と衝突する確率が p ( x ) = Σ t exp( −Σ t x ) なのではないことに注意。 x ~ x + dx で初めて媒質と衝突する確率が exp( −Σ t x )Σ t dx で与えられるのである。
また、④から 0 ~ X で衝突する確率は、X∫ Σ exp( −Σ x )dx = 1 − exp( −Σ X )
0

t

t

(4.8)

t

である。中性子は距離 X までに衝突するか、しないかどちらかであることから、この結果
は 1 から衝突しない確率である exp( −Σ t X ) を差し引いた結果になっている。□■例 8：平均値μ、標準偏差σの正規分布に対する確率密度関数は、次式で与えられる。

p( x ) =( x − µ )2
1
exp{−
}
2σ 2
2π σ

(4.9)

全ての変数区間で積分すると( x − µ )2
1∫− ∞ 2π σ exp{− 2σ 2 }dx = 1
例えば µ = 0, σ = 1 の標準正規分布の場合、∞P ( −1 < x < 1) = ∫1−1

P ( −2 < x < 2 ) = ∫2−2(4.10)x21
exp{− }dx = 0.6826
22π(4.11)x21
exp{− }dx = 0.954422π(4.12)x2
1(4.13)
exp{− }dx = 0.9974− 3 2π2(4.11), (4.12), (4.13)は正規分布の場合、µ ± σ に全体の約 68%、µ ± 2σ に約 95%、µ ± 3σ にP( −3 < x < 3) = ∫3

約 99.7%が含まれていることを示している。□
4.3 累積分布関数
ある確率分布関数 p (x ) が与えられたとき、 x < α である確率は、
α

F (α ) = P( −∞ < x < α ) = ∫ p( x )dx(4.14)−∞

13

で与えられる。これを累積分布関数という。例 7 で中性子が位置 0～X で衝突する確率を計算したが、これは累積分布関数になっている。
■例 8：偏差値は平均値が 50，
標準偏差が 10 の正規分布を仮定した場合の成績分布である。
偏差値の累積分布関数は、F (α ) = ∫α−∞( x − 50) 2
1
exp{−
}dx
2 × 102
2π × 10(4.15)

で与えられる。従って、偏差値 50 以下である確率はF (50) = ∫50−∞( x − 50) 2
1
exp{−}dx = 0.5
2 × 102
2π × 10(4.16)

偏差値 60 以下である確率はF (60) = ∫60−∞( x − 50) 2
1
exp{−}dx = 0.8413
2 × 102
2π × 10(4.17)

偏差値 90 以下である確率はF (90) = ∫90−∞( x − 50) 2
1
exp{−}dx = 0.9999683
2 × 102
2π × 10(4.18)である。(4.18)式から、偏差値 90 以上の人は(あくまでも成績が完全な正規分布になっている場合であるが)10 万人中約 3 人である。□

5.期待値、分散、標準偏差5.1 期待値
1～10 の数字が書かれたカードが 10 枚、裏返しておかれている。このカードを一枚引き、1 が出れば 1000 円もらい、2-10 が出れば 100 円払う。さて、あなたはこの賭け事に乗るか？
この問題を考えるとき、重要なのは、「どれぐらい得すると期待できるか」である。これは期待値を用いることで評価できる。
ある離散的なパラメータ X i に対して確率関数 p ( X i ) が与えられているとする。この場合、

X i の期待値はE[ X ] = ∑ X i p( X i )

(5.1)

i

で与えられる。これは、平均値である。
■例 9：先ほどの賭の期待値は1000 ×1
9− 100 × = 10 円
10
10

14

であり、平均的に得する賭になっている 8。□■例 10：さいころの目の期待値は

1×

1
1
1
1
1
1 7
+ 2 × + 3× + 4 × + 5× + 6 × =
6
6
6
6
6
6 2

である。□
ある連続的なパラメータ X に対して確率関数 p ( X ) が与えられているとする。この場合、X の期待値は∞E[ X ] = ∫ Xp( X )dX

(5.2)

−∞

で与えられる。
■例 11：巨視的断面積Σt の物質に中性子が入射したときに衝突せずに飛行する距離の期待値は確率密度関数が(4.6)式で与えられるため、∞

1∫ xΣ exp( −Σ x )dx = Σ
0

t

t

(5.3)

t

である。これは中性子が衝突なしに飛行する距離の平均を示したものであり、平均自由行程と呼ばれる。□
5.2 分散■例 12： 1～10 の数字が一つずつ書かれた 10 枚のカードが裏返しておかれている。この
カードを一枚引くが、二種類の賭ができる。
①1～5 が出れば 100 円もらい、6～10 が出れば 50 円払う。
②1～5 が出れば 1,000,001,000 円もらい、6～10 が出れば 1,000,000,500 円払う。
さて、あなたはどちらの賭に参加するか？①、②の期待値は①： 100 ×

55− 50 × = 25 円
10
10②： 1,000,001,000 ×55− 1,000,000,500 × = 250 円
1010であり、②の方が明らかに有利である。では、②の賭をやるだろうか？一般的には、②に8

少年マガジンに連載中の「賭博覇王伝 零」においては、これでは失格である。なぜならば、「1～10 の数字が書かれたカードが 10 枚、裏返しておかれている」と書かれているだけ
であり、1～10 の数字が書かれているカードがそれぞれ 1 枚ずつあるとは限らないからである。普通の人にとってはこの手の話は、予測がつかないために確かに「不確かさ」である。
一方、「賭博覇王伝 零」の主人公の宇海零にとっては、この手の話は「お見通し」であり、従って「不確かさ」ではなく、考慮すべき「系統誤差」なのかも知れない。
15

は手を出さないのではないかと思うが・・。□
例 12 は、ばらつきのあるデータを扱う場合、期待値のみでは情報として不十分であることを示している。すなわち、「どの程度データがばらついているか」を示す指標が必要である。分散は、データのばらつきを示す指標である。
ある離散的なパラメータ X = ( X 1 , X 2 ,..) に対して確率関数 p ( X i ) が与えられているとする。この場合、 X i の分散はvar[ X ] = E[( X − µ ) 2 ] = ∑ ( X i − µ ) 2 p( X i )

(5.4)

i

ここで、 µ = E [ X ] である。
■例 13：例 12 の①、②の賭の分散は
2① (100 − 25) ×

5
5+ ( −50 − 25) 2 × = 5625
10102② (1,000,001,000 − 250) ×55+ ( −1,000,000,500 − 250) 2 × ≈ 1018
1010であり、データのばらつきが全く異なる。つまり、②はばらつきが大きく「危険な賭」である。□■例 14：あるデータ列{1, 2, 3, 4, 5, 6}の分散は、このデータ列の平均が 3.5 であることから、(1 − 3.5) 2 + ( 2 − 3.5) 2 + (3 − 3.5) 2 + ( 4 − 3.5) 2 + (5 − 3.5) 2 + (6 − 3.5) 2≅ 2.9176である。□
ある連続的なパラメータ X に対して確率関数 p ( X ) が与えられているとする。この場合、X の分散は∞var[ X ] = E[( X − µ ) 2 ] = ∫ ( X − µ ) 2 p( X )dX
−∞

(5.5)

で与えられる。
■例 15：一様な確率密度関数がある。

1p( x ) = , − 1 < x < 12p( x ) = 0, otherwise(5.6)

この分散は、

16

∞

1

−∞

−1

var[ x ] = ∫ x 2 p( x )dx = ∫1 2
1
x dx =
2
3

(5.7)

である。□
5.3 標準偏差
標準偏差は、データのばらつきを表す指標の一つであり、分散の平方根として与えられる。

σ = var[ X ](5.8)

例 8 で示したように、データが正規分布に従ってばらついている場合、全データの約 68%が µ ± σ に、約 95%が µ ± 2σ に、約 99.7%が µ ± 3σ に含まれる。

5.4 モーメント
平均値と分散は、(5.2)式および(5.5)式からそれぞれ E [ X ] 、 E [( X − µ ) 2 ] で計算できる。
これらの形は、いずれも∞E[( X − c ) k ] = ∫ ( X − c ) k p( X )dX

(5.9)

−∞

と一般化できる。 E [( X − c ) k ] を c まわりの k 次モーメントと呼ぶ。1 次から 4 次までのモーメントは、以下のように呼ばれる。
平均：原点周りの 1 次モーメント
分散：平均周りの 2 次モーメント
歪度：平均周りの 3 次モーメントであり、分布の偏りを示す指標である。負であれば左にゆがんだ分布、ゼロであれば左右対称、正であれば右にゆがんだ分布となる。
尖度：平均周りの 4 次モーメントであり、値が大きいほど正規分布に比べて裾が広がっている(裾が重い)。

6.母集団の推定
6.1 母集団とは
我々が白鳥を「白い」と認識している。これは、すべての白鳥が白いと確認した結果ではなく、これまでに見た白鳥が白かったことから、「残りの白鳥も白いだろう」と推定した結果である。このように一部のサンプルから、残りの集団の性質を推定することを母集団の推定と呼ぶ 9。
9「ブラック・スワン―不確実性とリスクの本質」、ナシーム・ニコラス・タレブ著は、母集団推定の失敗について書かれた好著である。
17

母集団とは、調べたい要素すべてからなる集団であるが、一般的に母集団の要素の数は非常に多く、全てを調べることはできない。従って、限られた調査結果から母集団を推定する必要が生じる。母集団推定の例としては、世論調査や工業製品の抜き取り調査があげられる。
6.2 平均値の不確かさ(標準偏差)
モンテカルロコードを用いてある炉心体系の計算を初期乱数を変えて n 回行った。この結果を用いて、実効増倍率の平均値と標準偏差を以下の式から計算した。

Xi
i =1 nnµn = ∑(6.1)

σ n = var[ X i ] =( X i − µ )2
∑
n
i =1
n

(6.2)

この平均値と標準偏差は、n を非常に増やした場合に比べ、どの程度正しいだろうか？
まず、平均値について考えてみよう。実効増倍率が正規分布に従っていると仮定すると、各計算の無限増倍率は

Xi = µ ± σ(6.3)と書くことができる 10。n 回の計算の実効増倍率の平均値は

Xi
i =1 nnµn = ∑(6.4)

である。各計算の実効増倍率のばらつきが独立であるとすると、平均値のばらつき(平均値の標準偏差)を以下のように合成することができる 11。

σn =

2

σ
σ
σ 
  = n 2 =∑nni =1  n 
n

2

(6.5)

以上のことから、n 回の計算を抜き出して平均値を求めた場合、この平均値の不確かさ(平均値の標準偏差)は、

σ
n

となる。

■例 16：モンテカルロコードによる計算で得られる実効増倍率の真の標準偏差が 0.02 であったとする 12。この場合、n 回の計算で得られた実効増倍率の平均値に関する標準偏差は以10± を普通と異なった意味で使っていることに注意。 xi = µ + σ , µ − σ ではない。 xi が µを中心として標準偏差 σ でばらついているという意味である。
11

各項が

Xi
X
µ ±σ
σ
であるため、各項のばらつきを i =より と見積もることができる。
nnnn

12n 回の計算により得られた標準偏差でないことに注意。一般的には、母集団の標準偏差は18

下の通りである。
n=1 の場合：

0.02= 0.021

n=4 の場合：

0.02= 0.014n=9 の場合：

0.02≅ 0.006667
9n=100 の場合：

0.02
= 0.002100

計算回数が増えると、平均値の予測精度が上がっていくが、これは直感的にもそうだといえるだろう。□
6.3 標準偏差の推定
では次に、6.2 節で扱った問題において、n 回の計算結果から評価した標準偏差がどの程度正しいか検討してみよう。標準偏差の自乗である分散は、「それぞれの値と平均値との距

1 n∑ ( X i − µ )2 で計算すべきものである。
n i =1離の二乗を平均したもの」
を表している。従って、

ここで、µ は、モンテカルロ計算で得られた実効増倍率の母集団の平均値である。つまり、
モンテカルロ計算を非常に数多く実行したときに得られる平均値となる。しかし、今回は n回の計算のみを行っているので、 µ の値は分からない。そこで、次善の策として実施した n
回 の 計 算 で 得 ら れ た 実 効 増 倍 率の 平 均 値 µn を 使 わ ざ る を 得 な い。 つ ま り 、 分 散 を

1 n∑ ( X i − µn )2 で評価せざるを得ない。
n i =1f (c) =

1 n1 n2と言う関数は、
(
X
−
c
)
c
=
µ
=
∑ i
∑ X i の時に極小値を取る。しかし、
nn i =1n i =1µ ≠ µn であることから(つまり、n 回の実効増倍率の平均値と母集団の平均値が一致しない
ことから)、

1 n
1 n( X i − µn ) 2 は、 ∑ ( X i − µ ) 2 に比べて必ず小さくなる。つまり、n 回の∑
n i =1n i =1

計算で得られた実効増倍率から評価した分散は、必ず小さめに計算されると言うことである。
では、どの程度小さめに予測されるのだろうか。逆に、n 回の計算の実効増倍率の標準偏差から、実効増倍率の母集団の標準偏差は推定できるだろうか。
n 回の実効増倍率の分散は、以下のように計算できる。これを標本分散という。

( ) = n1 ∑ ( X − µ )標本分散： σ n

2

n

i =1

i

n

2

(6.6)

わからないが、ここでは簡単のため、母集団の標準偏差がわかっているものとしている。
19

一方、母集団の分散は、以下の式により推定可能である。これを不偏分散という。
不偏分散： (σ ) =2

1 n∑ ( X i − µn ) 2
n − 1 i =1

(6.7)

■例 17：例 16 で標本分散と不偏分散を比較する。実効増倍率を平均 1.01、標準偏差を 0.02、
サンプル数(計算回数)を n としてランダムに 100 セット、n 個の実効増倍率を「発生」させ、それぞれの標本分散と不偏分散を平均して求めた。結果を表 6.1 に示すが、標本分散は過小評価になっている一方、不偏分散にはこのような傾向は見られない。□
表 6.1 標本分散と不偏分散の比較
サンプル数5

標本分散より計算 不偏分散より計算した標準偏差した標準偏差0.01810.0202100.01900.020140

0.01990.0202以下では、不偏分散において、なぜ n でなく n-1 で割るのかを説明する。
(6.6)式のn∑ ( X − µ ) に着目する。これは、以下のように変形できる。
i

i =1

n

n

2

(σ ) = E[∑ ( X − µ ) ]
2

n

i =1

i

n

2

n

= E[∑{( X i − µ ) + ( µ − µn )}2 ]
i =1
n

n

i =1

i =1

n

= E[∑ ( X i − µ ) 2 + 2∑ ( X i − µ )( µ − µn ) + ∑ ( µn − µ ) 2 ]
i =1

n

n

i =1

i =1

n

n

i =1

i =1

= E [ ∑ ( X i − µ ) 2 − 2( µ n − µ ) ∑ ( X i − µ ) + n ( µ n − µ ) 2 ]= E[∑ ( X i − µ ) 2 − 2( µn − µ )( ∑ X i − nµ ) + n( µn − µ ) 2 ]n= E[∑ ( X i − µ ) 2 − 2( µn − µ )( nµn − nµ ) + n( µn − µ ) 2 ]i =1n= E [ ∑ ( X i − µ ) 2 − 2n ( µn − µ ) 2 + n ( µn − µ ) 2 ]i =1n= E [ ∑ ( X i − µ ) 2 − n ( µn − µ ) 2 ]

(6.8)

i =1

(6.8)式の右辺第一項は、以下のように計算できる。

20

nE[∑ ( X i − µ ) 2 ] = E[( X 1 − µ ) 2 ] + E[( X 2 − µ ) 2 ] + ... + E[( X n − µ ) 2 ]i =1= σ + σ + ... + σ = nσ
2

2

2

(6.9)

2

(6.9)式の右辺第二項は、以下のように計算できる。

E[n( µn − µ ) 2 ] = n × E[( µn − µ ) 2 ] = n ×

σ2
n

=σ2

(6.10)

従って、(6.8)式と(6.9)式を(6.7)式に代入し、
n

n

i =1

i =1

E [∑ ( X i − µn ) 2 ] = E [∑ ( X i − µ ) 2 − n( µn − µ ) 2 ] == nσ 2 − σ 2 = ( n − 1)σ 2(6.11)(6.11)式により

E[

1 n( X i − µn ) 2 ] = σ 2∑
n − 1 i =1

(6.12)1 n( X i − µn ) 2 の期待値が確かに母集団の分散 σ 2 になっていることが分かる。
∑
n − 1 i =1

であり、

すなわち、n 個のサンプルを抜き出して分散を計算する場合、(6.7)式を用いることで、より確からしく母集団の分散を推定できることが分かる。(6.7)式は、サンプル数 n が限定されることによる分散の過小評価を補正しているという観点から「不偏」分散と呼ばれる。なお、注意する必要があるのは、(6.7)式で得られる分散はあくまでも推定値ということである。つまり、不偏分散にも不確かさがあると言うことである 13。
6.4 信頼係数と信頼区間
解析コードを用いて核特性の解析を行った場合、真の値がどの範囲にあるのか、推定することは重要となる。例えば、臨界安全の分野において、使用している解析コードがどのような実効増倍率を算出するかを確認することは不可欠である。もちろん解析値には不確かさが含まれるため、ある確率である範囲に入る、という推定の仕方が望ましい。例えば、95%の確率で実効増倍率が 0.998～1.004 に入る、と言った推定である。この意味は、同じような推定を例えば 10000 回行った場合、約 9500 回は推定された範囲内に真の実効増倍率が存在するが、残り約 500 回については、この区間から外れている、ということである。このような推定を行う場合に用いる用語の定義を以下に示す。
信頼係数：どの程度の確率でその区間内に推定値が存在するかを示す係数。炉物理計算の分野では 95%がよく使用される。
信頼区間：信頼係数で与えられた確率で真の値が存在すると推定された範囲。
信頼限界：信頼区間の下限値と上限値。
不偏分散の分散は 2σ 4 /( n − 1) で与えられる。
(http://www.geocities.jp/ikuro_kotaro/koramu/335_mix.htm)

1321

■例 18：モンテカルロコードを用いてある炉心体系の計算を初期乱数を変えて 5 回行った。
得られた実効増倍率は 0.9986, 1.0282, 1.0211, 0.9804, 1.0022 であった。この炉心に対して、ヒストリー数を十分に増やした(=統計誤差を無視できる)計算を行った場合、得られる実効増倍率はどの程度になると推定できるだろうか。
実効増倍率の統計誤差は正規分布に従うため、実効増倍率の平均値も正規分布に従って変化する。実効増倍率の平均値は(3.1)式より、
0.9986 + 1.0282 + 1.0211 + 0.9804 + 1.0022= 1.0061
5

また標準偏差(不確かさ)は(3.2)式もしくは(6.7)式より(0.9986 − 1.0061) 2 +(1.0282 − 1.0061) 2 +(1.0211 − 1.0061) 2 +(0.9804 − 1.0061) 2 +(1.0022 − 1.0061) 2(5 − 1)= 0.0190

従って、平均値の標準偏差は(3.3)式より0.0190= 0.0085
5

では、統計誤差を含まない実効増倍率の計算値はが 1.0061 ± 1.96 × 0.0085 の範囲に入る
確率は 95%になる 14、と考えて良いだろうか。
ここで重要なのは、実効増倍率の平均値の標準偏差は 5 回の計算結果から推定された値であるということである。初期乱数を変化させて別に 5 回のモンテカルロ計算を行い、その結果から標準偏差を推定すると、0.0085 にはならないだろう、と言うことが想像できる。
つまり、平均値の標準偏差自体に不確かさが含まれていることになる。この不確かさを考慮するためには、平均値の標準偏差に不確かさがない場合に比べてより広い区間を考えなければならないことは明らかであろう。
実は、0.9986, 1.0282, 1.0211, 0.9804, 1.0022 という「実効増倍率」は、平均 1.01，標準偏差0.02 の正規分布から 5 回ランダムにサンプリングを行った結果である。従って、正確な(母集団の)標準偏差は 0.02 であり、これを「知っていた」とすると、5 回の平均値に対する標準偏差は、 0.02 / 5 = 0.0089 となる。従って、統計誤差を含まない実効増倍率の計算値が1.0061 ± 1.96 × 0.0089 の範囲に入る確率が 95%となる。
実際に、「実効増倍率」を 5 回発生させ、平均値と平均値の標準偏差を求め、(平均値±1.96
×0.0089)の中に真の平均値(1.01)が入っているかどうかを 10000 回繰り返して確認してみた。
この結果、(平均値±1.96×0.0089)の中に真の平均値が入っている回数は 9494 回であり、入っていなかった回数は 506 回であった。つまり、真の平均値が(平均値－1.96×平均値の標準偏差～平均値＋1.96×平均値の標準偏差)に入っていると予測した場合、95%の確率で正し
14

1.960

P( −1.960 < x < 1.960) = ∫−1.960x2
1exp{− }dx = 0.95 。(4.9)～(4.13)も参照。
22π22いといえるという結果になっている。
次に、5 回の実効増倍率の「計算結果」から平均値の標準偏差を(不偏分散を通じて)推定
し、(平均値±1.96×平均値の推定された標準偏差)の中に真の平均値が入っているかどうかを 10000 回繰り返し確認した。なお、10000 回の推定において、標準偏差の最大値は 0.0234、
最小値は 0.00029 となり、推定された標準偏差自体にかなりのばらつきがあることが確認できた。真の標準偏差は先ほど述べたように 0.0089 である。また、(平均値±1.96×平均値の推定された標準偏差)の中に真の平均値が入っているかどうかを確認した場合、 (平均値±
1.96×平均値の推定された標準偏差)の中に真の平均値が入っていた回数は 8824 回であり、入っていなかった回数が 1176 回であった。これは、真の平均値が存在する区間の見積もりを明らかに過小評価していることを示している。
以上の結果から、平均値の標準偏差を不偏分散により推定する場合には、この平均値の標準偏差の推定値自体に不確かさがあるため、推定区間をより大きく取らなければならないことが分かる。なお、参考までに 5 回の実効増倍率を 10000 セット発生させ、それぞれのセットについて、実効増倍率の平均値および平均値の標準偏差の分布を計算したものを
図 6.1 に示す。図 6.1 より、5 回の実効増倍率の平均値(真値=1.01)は 0.97～1.04 という広い範囲に分布し、実効増倍率の平均値の標準偏差(真値=0.0089)は 0.001～0.02 以上と大きくばらついている。
0.025

実効増倍率標準偏差

0.020

0.015

0.010

0.005

0.000
0.960

0.970

0.980

0.990

1.000

1.010

1.020

1.030

1.040

1.050

実効増倍率平均値

図 6.1 5 回の計算で得られた実効増倍率の平均値(横軸)と平均値の標準偏差(縦軸)の分布
ここで扱っている問題は、Student の t 分布と呼ばれる確率分布により考慮することができる。
t 分布を用いてサンプル数 n=5 に対する 95%の信頼係数に対する信頼区間を求めると、(平均値±2.776×平均値の推定された標準偏差)で与えられる。これは、先ほどの(平均値±
1.960×平均値の推定された標準偏差)に比べて明らかに広い範囲である。この区間内に真の平均値が入っているかどうかを上記と同様に 10000 回確認すると、入っていた回数は 9522
回、入っていなかった回数は 478 回であった。すなわち、95%の信頼係数に見合った適正な信頼区間を与えているといえる。□

23

ここで、Student 15の t 分布について簡単に説明しておく。パラメータ X 1 , X 2 ,..., X n が平均

値 µ 、標準偏差 σ の正規分布に従うとする。まず、(6.13)式で定義されるパラメータ α を考
える。

α=

X −µ
σ/ n

X=

X 1 + X 2 + ... + X n：平均値n(6.13)

σ / n ： X の真の標準偏差( X 1 , X 2 ,..., X n から不偏分散により推定したものではなく、母集団の標準偏差から計算したもの)
αは、(6.5)式を考慮すると、平均値 0、標準偏差 1 の正規分布になる。これより、信頼係数95% を 満 足 す る 信 頼 区 間 は − 1.96 <X − 1.96σ< µ < X + 1.96

nは、 X − 1.96

σ
n

σ
n

X −µ
< 1.96 、 す な わ ち 、σ/ nとなる。つまり、母集団の標準偏差が分かっている場合に< µ < X + 1.96σnにより、信頼係数 95%における X の信頼区間を推定することができる。
一方、母集団の標準偏差は一般には明らかでないため、不偏分散を通じて推定することとなる。この場合、

α'=

X −µ
U/ n

(6.14)

ここで、X=

X 1 + X 2 + ... + X n
：平均値
n
n

U=

∑( X − µ)
i =1

1

n −1

2

：不偏分散を通じて推定した母集団の標準偏差U：平均値の標準偏差の推定値
nStudent さん、とは妙な名前と思われるが、これはペンネームとのこと。William Sealy
Gosset さんがギネスビール社に勤めていたとき、ビールの品質管理の過程で t 分布についての知見を得たが、ギネスビール社は当時このような企業秘密を論文として公表することを許していなかった。そこで、Gosset さんは、Student というペンネームで t 分布について(会社にばれない形で)発表した。(Wikipedia より抜粋) ギネスを飲むときには、t 分布を思いだそう。
15

24

を対象とすることになるが、 α ' は、正規分布でなく、t(n-1)分布に従うことが知られている
16。n=5 の場合、信頼係数 95%を満足する信頼区間は

X − 2.776− 2.776 <X −µ< 2.776 、すなわち、
U/ 5

U
U
となる。2.776 は、t 分布表において、自由度= n-1 =5-1
< µ < X + 2.77655＝4、信頼係数 95%として与えられる値である。以下に代表的な α ' の値を示す。
n-1=4，

信頼係数 95%：2.776

信頼係数 98%：3.747

n-1=5，

信頼係数 95%：2.571

信頼係数 98%：3.365

n-1=10，信頼係数 95%：2.228

信頼係数 98%：2.764

n-1=50，信頼係数 95%：2.009

信頼係数 98%：2.403

n-1=240，

信頼係数 95%：1.970

信頼係数 98%：2.342

n-1=∞,

信頼係数 95%：1.960

信頼係数 98%：2.326

これより、n が多くなるほど、係数が小さくなり、n=∞の時には、 α ' は平均 0、標準偏差 1の正規分布の α と一致することが分かる。
6.5 推定臨界下限値増倍率
例 18 では、核計算コードの実効増倍率の平均値について、信頼区間の推定を行った。この例のように、核計算コードの平均誤差を知ることも重要である。一方、臨界安全の分野
では、「核計算コードで計算した実効増倍率がいくつ以下であれば未臨界と判断して良いか」と言う問いに応えなければならない。「未臨界と判断できる実効増倍率の計算値」を推
定臨界下限増倍率という。
仮に、母集団の平均値と分散が十分な精度でわかっていれば、推定臨界下限増倍率を求めることは容易である。例えば、母集団の平均がμ、分散がσ2 であれば、推定臨界下限増
倍率を µ − 2σ と設定すれば、臨界になる確率は、 (100%-95.44%)/2=2.28%以下になる。しかし、実際には、平均値と分散は標本(限られた数の計算結果)を通じて推定した値でしかな
く、それぞれに不確かさが含まれる。このような状況の下で、推定臨界下限増倍率を求めるのは簡単なことではない 17。推定臨界下限増倍率の評価方法については、別途 Appendixにまとめたので、そちらを参照して欲しい。
6.6 注意事項
限られた測定データから、母集団の確率分布を推定することは、もともと困難な問題である。炉心解析コードの不確かさの評価においては、「誤差が正規分布する」という仮定を

16

詳細については、例えば東京大学教養学部統計学教室編、「統計学入門」、東大出版会、
p.201 などを参照のこと。
17 具体的な計算手順については、例えば奥野浩、内藤俶孝、「臨界安全解析コードシステ
ム JACS の計算誤差評価」、JAERI-M 87-057 参照。
25

置く場合が多いが、対象とする解析体系を注意深く選ばないとこのような前提は成り立たない。例えば、熱中性子体系の UO2 燃料に対する実効増倍率の誤差分布と、高速中性子体
系の MOX 燃料に対する実効増倍率の誤差分布では、誤差の平均値も標準偏差も異なることが予想できる 18。従って、これらを同じグループとすると、正規分布を前提とした信頼区間
推定において誤った結果を与えることになる。
図 6.2 は、
連続エネルギーモンテカルロコード MVP と JENDL4.0 の組み合わせにより種々の臨界実験の実効増倍率を解析した結果である 19。実効増倍率の実験値はいずれも 1.0 であり、MVP の統計誤差は 0.02%Δk/k 程度と十分に小さい。
MVP/JENDL-4.0によるICSBEP解析結果の頻度分布
60
頻度

50

頻度

40

30

20

10

次の級

1.02301.02201.02101.02001.01901.01801.01701.01601.01501.01401.01301.01201.01101.01001.00901.00801.00701.00601.00501.00401.00301.00201.00101.00000.99900.99800.99700.99600.99500.99400.99300.99200.99100.99000.98900.98800.98700.98600.98500.98400.98300.98200.98100.98000.97900.97800.97700.97600.97500.97400.97300.97200

データ区間

図 6.2 MVP+JENDL4.0 による臨界実験解析結果
図 6.2 を見ると、MVP+JENDL4.0 の解析結果はおおむね正規分布をしているように見えるが、燃料の種類別に分けると、図 6.3 のようになる。図 6.3 を見ると、例えば U233 体系
は系統的に実効増倍率を過小評価しており、また、その分布も正規分布からは外れているように見える。従って、図 6.2 のデータによって MVP+JENDL4.0 の無限増倍率の信頼区間を評価した場合、U233 体系に対しては、適切な結果を与えないことが予想される。誤差が正規分布をなしていると仮定できるのは、あくまでも核特性が類似の体系のみであることに注意が必要である。

18

例えば、
「臨界安全評価関連コードの検証解析に関する報告書」、JNES/SAE09-031，原子力安全基盤機構、平成 21 年に様々な体系における実効増倍率計算誤差の頻度分布が示されている。
19 奥村啓介、長家康展、「JENDL-4.0 に基づく連続エネルギーモンテカルロコード MVP
用の中性子断面積ライブラリーの作成と ICSBEP ハンドブックの臨界性ベンチマーク解析
への適用」、JAEA-Data/Code 2011/010 の結果を名古屋大学の遠藤知弘助教が整理したもの。
26

MVP/JENDL-4.0によるICSBEP解析結果の頻度分布
60
U233: 64cases
PU: 241cases

50

MIX: 88cases
LEU+IEU+HEU: 255cases

頻度

40

30

20

10

1.02301.02201.02101.02001.01901.01801.01701.01601.01501.01401.01301.01201.01101.01001.00901.00801.00701.00601.00501.00401.00301.00201.00101.00000.99900.99800.99700.99600.99500.99400.99300.99200.99100.99000.98900.98800.98700.98600.98300.98500.98400.98200.98100.98000.97900.97800.97700.97600.97500.97400.97300.97200

データ区間

図 6.3 MVP+JENDL4.0 による臨界実験解析結果(燃料の種類別に分類)

7.共分散、相関係数、分散共分散行列
ここまでは、それぞれの(一系列の)データ列について平均や標準偏差の計算を行ってきたが、一般の不確かさ評価では、多変数を対象にして解析を行う。この場合、それぞれの変数の分布に加え、相互にどのような相関があるかということが不確かさの評価上重要になる。不確かさの合成については、次章で説明することとし、ここでは、変数間の相関を表す指標である共分散と相関係数を説明する。
二つのデータ列 X = ( X 1 , X 2 ,...) , Y = (Y1 , Y2 ,...) がある。これらのデータの共分散は

1 n∑ ( X i − X )(Yi − Y )n i =1

cov[ X ,Y ] = E[( X − X )(Y − Y )] =

(7.1)

として与えられる。また、

α=E[( X − X )(Y − Y )]2=

2

E[( X − X ) ] E[(Y − Y ) ]cov[ X ,Y ]

σ[X ]

2

σ [Y ]

2=

cov[ X ,Y ]σ [ X ]σ [Y ](7.2)

を相関係数とよぶ。 α が-1 に近くなるほど強い負の相関を表し、逆に 1 に近くなるほど強い正の相関を表す。また、0 の場合は無相関であることを示す。例えば、 yi = xi であると
すると、

α=

E[( X − X )( X − X )]E[( X − X ) 2 ] E[( X − X ) 2 ]

=1

(7.3)

Yi = − X i であるとすると、

α=E[( X − X )( − X + X )]2

2

E[( X − X ) ] E[( − X + X ) ]=−E[( X − X )( X − X )]E[( X − X ) 2 ] E[( X − X ) 2 ]なお、定義から cov[ X , X ] = var[ X ] であることに注意。

27= −1 (7.4)

次に多くのデータ列 X 1 , X 2 ,..., を考える 20。これらのデータ列のうち、 X I と X J の共分
散は、

cov[ X I , X J ] = E[( X I − X I )( X J − X J )] =1 n
∑ ( X I − X I )( X J − X J )n i =1(7.5)

ここで

X I = E[ X I ] 、 X J = E[ X J ]である。これらの共分散を以下のように行列形式に並べたものが分散共分散行列である。

 cov[ X 1 , X 1 ] cov[ X 1 , X 2 ] cov[ X , X ] cov[ X , X ]
2
1
2
2
C=
..
..


cov[ X m , X 1 ] cov[ X m , X 2 ].. cov[ X 1 , X m ] .. cov[ X 2 , X m ] 

..
..


.. cov[ X m , X m ](7.6)

分散共分散行列の対角要素は、データ列 X I の分散(標準偏差の自乗)となり、非対角要素は、データ列 X I と X J の共分散となる。定義上、分散共分散行列は、対称行列である。また、

 cov[ X 1 , X 1 ] σ [ x ]σ [ x ]
1
1

X
X
cov[
,
2
1]

 σ [ x2 ]σ [ x1 ]..
 cov[ X m , X 1 ] σ [ xm ]σ [ x1 ]cov[ X 1 , X 2 ]σ [ x1 ]σ [ x2 ]cov[ X 2 , X 2 ]σ [ x2 ]σ [ x2 ]..
cov[ X m , X 1 ]σ [ xm ]σ [ x2 ]cov[ X 1 , X m ] 
σ [ x1 ]σ [ xm ] cov[ X 2 , X m ] ..
σ [ x2 ]σ [ xm ] ..
..
cov[ X m , X m ] ..
σ [ xm ]σ [ xm ] 

..

(7.7)

は相関行列と呼ばれ、各要素の相関係数を示したものとなる。
■例 19：2012 年 1 月 4 日～2012 年 6 月 27 日までの間の TOPIX(東京株価指数)、金価格、為替レート(円/ドル)の変動の相関を評価する。期間中のこれらの値の変化を図 7.1 に、相関関係を図 7.2 に示す。

20

X 1 = { X 1,1 , X 1, 2 ,..., X 1,i ,...}, X 2 = { X 2,1 , X 2, 2 ,..., X 2,i ,...},... である。
28

900

4800

850
800

4600

750

金価格

4400

700
4200

650

金価格
TOPIX
円/ドル(×10)

4000
3800

600
550

2012/6/20
2012/6/27

2012/5/23
2012/5/30
2012/6/6
2012/6/13

2012/5/2
2012/5/9
2012/5/16

2012/4/25

2012/4/4
2012/4/11
2012/4/18

2012/3/7
2012/3/14
2012/3/21
2012/3/28

2012/2/15
2012/2/22
2012/2/29

2012/1/25
2012/2/1
2012/2/8

2012/1/18

500

2012/1/4
2012/1/11

3600

TOPIXもしくは円/ドル(×10)

5000

図 7.1 TOPIX、金価格、為替レートの変動
900

900

850

850

800

800

4900
4800

750
700

4600

金価格

トピックス終値

トピックス終値

4700

750
700

4500
4400
4300
4200

650

650

4100
600
4000

600
4100

4200

4300

4400 4500
金価格(東京)

図 7.2

4600

4700

4800

4900

4000
75

77

79

81

83

85

75

77

79

円/ドル

81

83

85

円/ドル

TOPIX(東京株価指数)、金価格、為替レート(円/ドル)の相関関係

また、これら 3 つのパラメータの分散共分散行列および相関行列を表 7.1 に示す。□
表 7.1

TOPIX(東京株価指数)、金価格、為替レートの分散共分散行列および相関行列
分散共分散行列金

トピックス終値

金32311

トピックス終値

7854

2312

円/ドル

209

71.4

相関行列
円/ドル

金

トピックス終値

金1.00

トピックス終値

0.91

1.00

円/ドル

0.58

0.74

4.00

円/ドル

1.00

8.不確かさの伝搬
8.1 概要
ばらつきをそれぞれ有する二つのデータ列 X i , Yi がある。これらのデータ列の標準偏差は

σ x 、 σ y であるとする。 X i + Yi のばらつきはどのように評価できるだろうか。ここで、以下の三ケースについて考えてみよう。
① X i , Yi のばらつきが全く独立な場合② X i , Yi のばらつきに強い正の相関がある場合③ X i , Yi のばらつきに強い負の相関がある場合①、②、③どのケースについて、 X i + Yi のばらつきが大きくなるだろうか。②は、 X i の29

ばらつきが正側の時、 Yi のばらつきも正側になるため、ばらつきが足し合わされ大きくな
りそうである。最も極端なケースとして、 Yi = X i であれば、 X i + Yi の標準偏差は 2σ x にになる。一方、③は X i のばらつきが正側の時、 Yi のばらつきが負側になるため、ばらつき
が打ち消しあい、小さくなりそうである。この極端なケースは、 Yi = − X i であり、 X i + Yiの標準偏差は零になる。①は、②と③の間に来ると推定される。
①、②、③のケースについて、 X i + Yi の標準偏差は以下のように与えられる。これらの導出については後述する。
① σ x2 + σ y2

(8.1)

②σ x + σ y(8.2)③ σx −σ y(8.3)(8.1)式～(8.3)式を見ると、不確かさの大きさとしては、確かに②＞①＞③となっている。以下では、より一般的な形でこのような不確かさの伝搬について説明する。
8.2 多変数関数の誤差伝搬
ある多変数の関数 f ( X 1 , X 2 ,..., X m ) がある。 X 1 , X 2 ,..., X m の不確かさ ∆X 1 , ∆X 2 ,..., ∆X mに起因する f ( X 1 , X 2 ,..., X m ) の変化 ∆f は、(8.4)式のように表すことができる。

∆f =

m
∂f
∂f
∂f
∂f
∆X 1 +∆X 2 + ... +∆X m = ∑
∆X k
∂X 1
∂X 2
∂X m
k =1 ∂X k(8.4)(8.4)式の両辺を二乗する。

(∆f )2 = ∑∑ ∂f∂f∆X l ∆X k
l =1 k =1 ∂X l ∂X k
m

m

(8.5)

(8.5)式の両辺の平均(期待値)を計算する 21。

E[(∆f ) ] = E[( f − f ) ]
2

2

m m
∂f ∂f
∂f ∂f
∆X l ∆X k ] = ∑∑E[ ∆X l ∆X k ]l =1 k =1 ∂X l ∂X k
l =1 k =1 ∂X l ∂X kmm

= E[∑∑

(8.6)

∂f ∂f
E[( X l − X l )( X k − X k )]l =1 k =1 ∂X l ∂X k
m

m

= ∑∑

なお、 f = E [ f ( X 1 , X 2 ,..., X m )] 、 X l = E [ X l ] 、 X k = E [ X k ] である。
21

変数が一つの際の平均値の計算は(5.2)式であった。多変数関数の場合には
∞

∞

−∞

−∞ −∞

E[ f ( X 1 , X 2 ,..., X m )] = ∫ ...∫∞∫ f ( X , X ,..., X ) p( X , X ,..., X )dX dX ,..dX
1

2

となる。
30

m

1

2

m

1

2

m

(8.6)式より、(8.7)式が得られる。

∂f ∂fE[( X l − X l )( X k − X k )]l =1 k =1 ∂X l ∂X kmm

E[( f − f ) ] = ∑∑
2

(8.7)

E[( f − f ) ] は f ( X 1 , X 2 ,..., X m ) の分散を、E[( X l − X l )( X k − X k )] は X l と X k の共分散を2示している。従って、 f ( X 1 , X 2 ,..., X m ) の分散は、変数の共分散 E [( X l − X l )( X k − X k )] と
偏微分値

∂f ∂f
∂f
から計算できる。なお、は X l が変化したときの f の変化量を示し,∂X l ∂X k∂X lており、 f に対する X l の「効き方」を表しているので感度係数と呼ばれる 22。
8.3 二変数の場合の誤差伝搬では、(8.7)式を二変数の場合について詳しく見てみよう。二変数の場合、(8.7)式は、E [( ∆f ) 2 ] =

∂f ∂f
∂f ∂f
E [ ∆X 1∆X 1 ] +E [ ∆X 1 ∆X 2 ]∂X 1 ∂X 1∂X 1 ∂X 2

(8.8)

∂f ∂f
∂f ∂f
+
E [ ∆X 2 ∆X 1 ] +E [ ∆X 2 ∆X 2 ]∂X 2 ∂X 1∂X 2 ∂X 2

X 1 , X 2 の ば ら つ き に 全 く 相 関 が な い 場 合 ( 共 分 散 =0 の 場 合 ) 、E[ ∆X 1∆X 2 ] = E[ ∆X 2 ∆X 1 ] = 0 であることから、(8.8)式は(8.9)式の形に書くことができる。
E[( ∆f ) 2 ] =

∂f ∂f
∂f ∂f
E[ ∆X 1∆X 1 ] +E [ ∆X 2 ∆X 2 ]∂X 1 ∂X 1∂X 2 ∂X 2(8.9)

両辺の平方根を取って

σ ( f ( X 1 , X 2 )) = var( f ( X 1 , X 2 )) = E [( ∆f ) 2 ]
=

∂f ∂f
∂f ∂f
E [ ∆X 1∆X 1 ] +E [ ∆X 2 ∆X 2 ]∂X 1 ∂X 1∂X 2 ∂X 222

 ∂f 
 ∂f 
 cov[ X 2 , X 2 ] cov[ X 1 , X 1 ] + 
= 
∂
∂
X
X
 2
 1
2

(8.10)

2 ∂f 
 ∂f 
 var[ X 2 , X 2 ] var[ X 1 , X 1 ] + 
= 
∂
∂
X
X
 2
 1
2


  ∂f ∂f
σ ( X 1 )  + σ ( X 2 ) =   ∂X 2
 ∂X 1
22

2

X ∂f
∂f
=
f ∂X lf

炉物理分野の感度解析で出てくる感度係数は、相対感度係数であり、 l

∂X l
Xl

で表される。これは、 X l の相対的な変化量に対する f の相対的な変化量を示している。例
えば、 X l が 1%変化したときの f の変化量を%単位で示したものである。
31

であり、通常の誤差の合成になる。
■ 例 20 ： X 1 , X 2 の 不 確 か さ に よ る f1 ( X 1 , X 2 ) = X 1 + X 2 、 f 2 ( X 1 , X 2 ) = X 1 X 2 、

f 3 ( X 1 , X 2 ) = X 1 / X 2 の不確かさ(標準偏差)を計算する。なお、 X 1 , X 2 の不確かさに相関はないものとする。(8.10)式を用いると、

σ ( f1 ( X 1 , X 2 )) =(σ ( X 1 ) )2 + (σ ( X 2 ) )2

σ ( f 2 ( X 1 , X 2 )) =( X 2σ ( X 1 ) )2 + ( X 1σ ( X 2 ) )2

(8.11)

2= X1X 2 σ ( X1)   σ ( X 2 )  + 
 X1   X 2 (8.12)

2

2


 
 1
X1
σ ( f 3 ( X 1 , X 2 )) =  σ ( X 1 )  +  −σ ( X 2 ) 2
  (X 2 )
 X2

=

X1
X2

2

 σ ( X1)   σ ( X 2 )  + 
 X1   X 2 22(8.13)

以上。□
次に、 X 1 , X 2 の相関係数が α である場合を考える。(7.2)式より

α=E[ ∆X 1∆X 2 ]E[ ∆X 1∆X 1 ] E[ ∆X 2 ∆X 2 ](8.14)

従ってE[ ∆X 1∆X 2 ] = α E[ ∆X 1∆X 1 ] E[ ∆X 2 ∆X 2 ] = ασ ( X 1 )σ ( X 2 )(8.15)(8.8)式より

σ ( f ( X 1 , X 2 )) =

∂f ∂f
∂f ∂f
E [ ∆X 1∆X 1 ] +E [ ∆X 1∆X 2 ]∂X 1 ∂X 1∂X 1 ∂X 2
+

∂f ∂f
∂f ∂f
E [ ∆X 2 ∆X 1 ] +E [ ∆X 2 ∆X 2 ]∂X 2 ∂X 1∂X 2 ∂X 222  ∂f ∂f∂f ∂f
σ ( X 1 )  + σ ( X 2 )  + 2α
σ ( X 1 )σ ( X 2 )= ∂X 1 ∂X 2  ∂X 2
 ∂X 1

α = 1 の場合

32

(8.16)

22  ∂f
 ∂f∂f ∂f
σ ( f ( X 1 , X 2 )) = σ ( X 1 )  + σ ( X 2 )  + 2
σ ( X 1 )σ ( X 2 )∂X 1 ∂X 2  ∂X 2
 ∂X 1
2(8.17)


 ∂f
∂f
∂f
∂f
σ ( X1) +σ ( X 2 )  =σ ( X1) +σ (X2)
= 
∂X 1
∂X 2
∂X 2

 ∂X 1

α = −1 の場合
22  ∂f
 ∂f∂f ∂f
σ ( f ( X 1 , X 2 )) = σ ( X 1 )  + σ ( X 2 )  − 2
σ ( X 1 )σ ( X 2 )∂X 1 ∂X 2  ∂X 2
 ∂X 1
2(8.18)


 ∂f
∂f
∂f
∂f
σ ( X1) −σ ( X 2 )  =σ ( X1) −σ (X2)
= 
∂X 1
∂X 2
∂X 2

 ∂X 1■例 21： f ( X 1 , X 2 ) = X 1 + X 2 であれば、∂f∂f== 1 であることから、本章の最初に∂X 1 ∂X 2

示した①、②、③のそれぞれのケースが正しいことが分かる 23。□

8.4 不確かさ伝搬の Sandwich formula
2
次に、核特性の感度解析で用いる不確かさ伝搬の式を導いておく。(8.7)式を f で割るこ

とにより、
2

m m
 ∆f 
X X ∂f ∂f∆X ∆X kE[  ] = ∑∑ l 2 k
E[ l]f ∂X l ∂X k
Xl Xk
l =1 k =1
 f 
m m∂f / f∂f / f∆X ∆X k
E[ l
]
= ∑∑
Xl Xk
l =1 k =1 ∂X l / X l ∂X k / X kmm= ∑∑ S f , x l S f , x k E[
l =1 k =1

ここで S f , X l =(8.19)∆X l ∆X k]
Xl Xk∂f / f∂f / fである。これらは、脚注 15 で述べたように相対,S f ,Xk =∂X l / X l∂X k / X k

感度係数である。(8.19)式を行列形式で書くと

23

誤差をベクトルと考え、共分散をこれらのベクトルのなす角と考えるとわかりやすい。
共分散(相関係数)=0 は二つの誤差ベクトルが直角であり、二つのベクトルの大きさは
22x + y で合成される。共分散(相関係数)=1 は二つの誤差ベクトルが平行であり、二つのベクトルの大きさは x + y で合成される。
33

2

[

 ∆f 
E [  ] = S f , X 1
 f 

S f ,X2∆X 1 ∆X 2
 ∆X 1 ∆X 1
E
E
[
]
[
]

X1 X1
X1 X 2
 E [ ∆X 2 ∆X 1 ] E[ ∆X 2 ∆X 2 ]... S f , X m 
X 2 X1
X2 X2

..
..
 ∆X m ∆X 1
∆X m ∆X 2
] E[
]
 E[
X
X
X
X
1
2
m
m


]

∆X 1 ∆X m ]
X1 X m  S f , X1 ∆X ∆X m   S.. E[ 2]  f ,X2 
X 2 X m   ... 

..
..

∆X ∆X m   S f , X m .. E [ m]Xm Xm ..

E[

(8.20)である。これはC ff = t S f C xx S f(8.21)と表記することができる。 C ff は f の相対分散共分散行列(行列サイズ：1×1)、 C xx はX 1 , X 2 ,..., X m の相対共分散行列(行列サイズ：m×m)、 S f は f の相対感度行列(行列サイズm×1)である。(8.19)式は、共分散行列を感度行列で「挟み込む」形になっているため、炉物理の感度解析業界では"Sandwich Formula"などと呼ばれている。(8.20)式の平方根を取ること
で f の相対標準偏差を求めることができる。
次に考慮すべき関数が複数 M 個ある場合を考える。この場合、

 ∂f1 / f1
 ∂X / X
 1 1
 ∂f1 / f1
S f =  ∂X 2 / X 2...
 ∂f / f
 1 1
 ∂X m / X m∂f 2 / f 2∂X 1 / X 1∂f 2 / f 2∂X 2 / X 2∂f 2 / f 2∂X m / X m∂f M / f M ∂X 1 / X 1 ∂f M / f M ∂X 2 / X 2 

...
...
∂f M / f M ...
∂X m / X m ...

となり、S f のサイズは m×M、C xx のサイズは m×m、tS のサイズは M×m であるので、C ff は M×M の行列になる。 C ff の対角要素は、 X 1 , X 2 ,..., X m の不確かさに起因するそれぞれの関数の相対値の分散、非対角要素は断面積の不確定性に起因する二つの関数の相対値の共分散を示している。
■例 22： X 1 , X 2 ,..., X m が多群の微視的断面積、 f i が少数群の集合体平均巨視的断面積であるとする。この場合、 S f は、多群の微視的断面積に対する少数群集合体平均巨視的断面積の感度係数を表す行列となる。この感度係数行列を何らかの方法で求めることができれば、

34

(8.19)式により少数群集合体平均巨視的断面積の分散と、共分散を求めることができる。この分散は、例えば 1 群の吸収断面積、2 群の核分裂断面積の不確かさを示すものである。また共分散は、1 群の吸収断面積、2 群の核分裂断面積がどの程度の相関を持つかを表すことになる。□■例 23： x1 , x 2 ,..., x m が少数群の集合体平均断面積、 f i が核特性値(例：実効増倍率、出力分布、反応度係数、・・・)であるとすると、(8.21)を用いることで少数群集合体平均断面積の不確かさに起因して炉心核特性値にどの程度の不確かさを生じるか、また炉心特性値間にどの程度の相関が生じるのかを計算することができる。□

9.まとめ
不確かさの評価技術は、解析技術の進歩に伴ってより重要になる。しかし、不確かさ評
価で必要となる計算は多くの場合複雑な統計処理の理論に基づいており、これらを理解す
ることはそれなりにハードルが高い。そこで、本稿では、この不確かさ評価を実施する上で必要となる基礎的な事項をまとめ、不確かさ評価に取り組む上での手助けになる参考資料とすることを目指した。

謝辞
推定臨界下限増倍率の評価(Appendix)については、名古屋大学の遠藤知弘氏にまとめて頂いたものである。また、本稿の草稿について、大阪大学の北田孝典氏、北海道大学の千葉豪氏、名古屋大学の遠藤知弘氏に有益なコメントを頂いた。ここに記して感謝の意を表したい。

35

Appendix. 推定臨界下限増倍率の算出法
本付録では、推定臨界下限増倍率 k sub の算出方法について補足説明を行う。ここでは n個の臨界実験を解析して、実効増倍率 k eff の標本として ki (i=1,…n)を得たとして話を進める。
2

得られた n 個の k eff の解析結果から、k eff の標本平均 k および不偏分散 s を以下のように推定することができる。

k=

1 n
∑ ki
n i =1

(A.1)

2
1 n
(
ki − k )∑n − 1 i =1

(A.2)

s2 =

ここで k eff の母集団は平均を µ 、標準偏差 σ の正規分布に従うと仮定する。

■例題：あるモンテカルロコードで 10 個の臨界実験を解析したところ、
1.0412, 1.0236, 1.0198, 1.0297, 1.0087, 0.9741, 1.0195, 1.0200, 1.0243, 0.9818という値が得られた。このモンテカルロコードの誤差が正規分布していると仮定できる場
合、推定臨界下限増倍率を求めよ。(言い換えると、実効増倍率がいくら以下であれば未臨界と言えるかを明らかにせよ、ということである。)
なお、本計算の平均値 k と標準偏差 s は以下のように求められる。
k=

1.0412 + 1.0236 + 1.0198 + 1.0297 + 1.0087 + 0.9741 + 1.0195 + 1.0200 + 1.0243 + 0.9818= 1.014310(1.0412 − 1.0143 ) 2 +(1.0236 − 1.0143 ) 2 +(1.0198 − 1.0143 ) 2 +(1.0297 − 1.0143 ) 2 +(1.0087 − 1.0143 ) 2s=+ (0.9741 − 1.0143 ) 2 +(1.0195 − 1.0143 ) 2 +(1.0200 − 1.0143 ) 2 +(1.0243 − 1.0143 ) 2 +(0.9818 − 1.0143 ) 2(10 − 1)

= 0.0209

A.1 必要となる統計学知識
推定臨界下限増倍率を算出するに当たって、必要となる統計学の知識を以下で記述する。
(1) 標本平均 k は、平均 µ 、標準偏差 σ

n の正規分布に従う。

2(2) 不偏分散 s については、

z = (n − 1)

s2

σ2

n

=∑
i =1

(k − k )

2

i

(A.3)

σ2
36

という量が自由度 (n − 1) のカイ二乗分布に従う。カイ二乗分布の性質より、z の平均は

n − 1 、分散は 2 (n − 1) となる。従って不偏分散 s 2 について考えた場合には、 s 2 の平均は
σ 2 、分散 2 σ 4 (n − 1) となる。n が十分に大きくなるにつれて、 s 2 は平均 σ 2 、標準偏差2 (n − 1) σ 2 の正規分布に近づくようになる。

(3) 変数 X として平均 λ、分散 1 の正規分布に従う量を考える。また、変数 Y として自由度

ν のカイ二乗分布に従う量を考える。ここで X と Y が統計的に独立である場合、以下で定義される量 T は、自由度ν 、非心パラメータ λ の非心 t 分布に従う。
T=

X
Yν

(A.4)

A.2 実効増倍率の母集団の平均値と標準偏差がわかっている場合
まず、最も簡単なケースとして、実効増倍率の母集団の平均 µ と標準偏差 σ がわかっている場合を考える。例題として挙げた 10 個の実効増倍率の結果は、実は平均 1.01、標準偏差 0.02 の正規分布に従って、10 回ランダムにサンプリングを行った結果である。すなわち

k eff の母集団は µ = 1.01 、標準偏差 σ = 0.02 の正規分布に従っている。この母集団の平均と標準偏差を用いると、例えば 1.01 − 2 × 0.02 < k eff < 1.01 + 2 × 0.02 になる確率は 95.44%
で あ る こ と が 言 え る 。 こ れ は 、 k eff < 1.01 − 2 × 0.02 に な る 確 率 ( 下 側 累 積 確 率 ) が2.28%(=(100%-95.44%)/2)であることと等価である。従って、推定臨界下限増倍率 k sub として、
例えば 0.97(=平均値－2×標準偏差=1.01－2×0.02)にした場合、 k eff が 0.97 未満であれば、臨界にならない確率は 100%-2.28%=97.72%となる。逆に、臨界になる確率は 2.28%である(図A.1)。これは、 k eff が 0.97 未満となる解析値が 10000 ケースあった場合、そのうちの 228 ケースが臨界になることを示している。判断基準を 0.95(=1.01－3×0.02)とした場合、
つまり k eff
が 0.95 未満であれば、臨界にならない確率は先と同様の手順により 99.87%になる。以上で議論したように、臨界になる確率を零にはできないが、推定臨界下限増倍率 k sub を小さく設定することにより、限りなく小さくできることを示している。

37

keff
頻度分布

2σ
この割合
が2.28%

図 A.1 母集団の平均と標準偏差がわかっている場合の k sub ( µ , σ ) 評価例
以上で述べた手順をまとめると、母集団の平均 µ と標準偏差 σ がわかっている場合については k sub ( µ , σ ) を以下のような手順で設定することができる。
①臨界になる確率 p を設定する。
②実効増倍率 k eff が平均 µ 、標準偏差 σ の正規分布に従うとして、正規分布の下側累積確率 p に対応するパーセント点 − α p を求める。すなわち k eff の分布を k eff ≤ µ − α p σ まで積分した累積確率が p に対応する。
③推定臨界下限増倍率を k sub ( µ , σ ) = µ − α p σ と設定する。このように設定することで、

k eff の解析結果が k sub ( µ , σ ) 未満であれば、臨界にならない確率が (1 − p ) となる。
以上の手順に基づいて例題の推定臨界下限増倍率 k sub ( µ , σ ) を算出してみると、臨界になる確率を p = 2.5% と設定した場合には、それに対応するパーセント点は α p = 1.96 となる
ため、 k sub ( µ , σ ) = 1.01 − 1.96 × 0.02 = 0.9708 となる。

A.3 実効増倍率の母集団の平均値と標準偏差が不明であるが、これらが十分に精度良く推定できている場合
しかし、実際には母集団の平均 µ と標準偏差 σ は未知であり、これらは計算結果(標本)2の平均 k と不偏分散 s により推定される。ここで重要なことは、標本を用いた平均 k も不2

偏分散 s もいずれも不確かさを含んでいると言うことである。つまり、標本を取り直すと(別の体系を計算してみると)、平均と不偏分散としては、別の値が得られるということである。ただし、このような状況であったとしても、標本数 n が十分に多い場合には、平均と38

不偏分散の不確かさを十分小さくすることができる。これは A.1 節で述べたように、標本数n → ∞ の極限で k と s 2 の分散(標準偏差)が零に漸近することからも理解できる。従って、標本数が十分に大きく、平均と不偏分散の不確かさが十分に小さければ、図 A.2 に示すように、図 A.1 と(厳密には同等でないが)ほぼ同じようなやり方で推定臨界下限増倍率を求める
ことができそうである。

k − 2s
確率分布

keff
頻度分布

2σ
この割合
が2.28%

図 A.2 母集団の平均と標準偏差が未知、標本数が十分に大きい場合の k sub 評価例
従って、標本数が十分に大きく k と s の推定結果が母集団の平均 µ と σ に限りなく近い
22

場合については、推定臨界下限増倍率 k sub を以下のような手順で設定できそうである。
2①n 個の実効増倍率 k eff の解析結果から、平均 k と不偏分散 s を求める。
②臨界になる確率 p を設定し、実効増倍率 k eff が平均 µ 、標準偏差 σ の正規分布に従うとして、正規分布の下側累積確率 p に対応するパーセント点 − α p を求める。
③推定臨界下限増倍率を k sub = k − α p s と設定する。

以上の手順に基づいて例題の推定臨界下限増倍率 k sub を算出してみると、 p = 2.5% と設定した場合には、 k sub = 1.0143 − 1.96 × 0.0209 = 0.9733 となる。この値は A.2 節で求められた、母集団の平均 µ と σ から推定された k sub ( µ , σ ) の値 0.9708 に比較的近いものの、若干大きな値となっている。こうして求められた k sub の値 0.9733 は、A.2 節で述べた µ と σ を用いた評価方法で考えた場合には、 k sub ( µ , σ ) = 0.9733 = 1.01 − 1.835 × 0.02 に対応しており、実質的には α p の値を小さめに設定してしまったことに対応する。すなわち、臨界になる確率が p = 2.5% となるように k sub を設定したつもりであったが、µ と σ のかわりに用いた k と s の不確かさによって、臨界になる確率は 3.3%とわずかに増えてしまっている。

39

例題で示された 10 個の k eff 解析結果の場合は、標本数 n が 10 個でも k sub を(たまたま)そ
れなりの値で設定することができた。参考までに、10 個で 1 セットの臨界実験の解析を別途 1000 ケース実施したことを想定してランダムサンプリングを行い、本節で述べた手法により各ケースで k sub を算出した結果を図 A.3 に示す。図 A.3 に示したように、全 1000 ケー
スのうち約 56%が k sub ( µ , σ ) = 0.9708 よりも大きく設定される結果となり、これら約 56%のケースについては、「臨界になる確率が p = 2.5% 超えることがないように k sub を設定できていない」こととなる。つまり、ここで示した手順が適用できるのは、あくまでも母集団の平均と標準偏差が「きわめてよく推定できている」場合のみであり、サンプル数が 10

推定臨界下限増倍率

程度のケースに使用することはできないことがわかる。

1.02
1.01
1.00
0.99
0.98
0.97
0.96
0.95
0.94
0.93
0

200

400

600

800

1000

ケ ース
図 A.3 k − α p s により設定した k sub の散布図(標本数 n = 10 , p = 2.5% )A.4 実効増倍率の母集団の平均値が不明であり、標準偏差のみがわかっている場合
A.3 節では µ と σ のかわりとして k と s を用いた k sub の設定方法について述べたが、現実的には、 k と s が µ と σ に等しいとみなせるほど十分な標本数を得ることはできない場合もある。この場合、後述の図 A.6 に示すように、(推定臨界下限増倍率の候補となり得る)「平均値－ α ×標準偏差( k − α s )」の値自体がある範囲でばらつくことになる。本節では、次のステップとして、実効増倍率の母集団の平均値が不明で( k のみ不確かさがあり)、母集団の標準偏差がわかっているという仮想的な場合を考えることとする(図 A.4)40

k −α σ

keff

確率分布

頻度分布

ασ
この割合
がγ%

この割合
がp%

図 A.4 母集団の平均値が不明であるが、標準偏差が分かっている場合

この場合、標本平均 k は平均が µ 、標準偏差が σn の正規分布に従ってばらつくことになる。ここで仮に、実効増倍率の平均値 k が真値 µ より高めに推定されたとする。例えば、1.01 + 2

0.02= 1.0226 という平均値が得られたとする(1.0226 より高い平均値 k が得ら10れる確率 γ は 2.28%である)。ここで、臨界になる確率 p を 2.28%となることを期待して、

α = 2 として推定臨界下限増倍率 k sub を k − α σ = 1.0226 − 2 × 0.02 = 0.9826 と推定してみる。こうして得られた値は、A.2 節で述べた µ − 2σ で求められる臨界下限増倍率 0.97(=1.01－2×0.02)より明らかに大きい。すなわち、 k sub として 0.9826 を設定してしまうと臨界になる確率は 2.28%より大きくなる(この場合、実質的には α p = 1.37 と設定したこととなり、臨界になる確率は 8.53%と大きくなる)。従って、 k のばらつきまで考慮して、臨界になる確率が 2.28%超えないようにするためには、 α の値を 2 よりも大きくとる必要がある。
A.2, A.3 のケースでは、目標になる実効増倍率 k − α σ が幅を持たない「線」だった。従って、どれぐらいの確率で臨界になるか、ということだけを考えていれば良かった。しかし、今回は、目標になる実効増倍率 k − α σ が「幅を持つ領域」になってしまっている。ということは、この「幅を持つ領域」の中で、どこに目標値(つまり、どれぐらいの確率で臨界になるかというしきい値の線)を置くか、ということも追加的に考えないといけないことになる。考えないとならないパラメータが二つに増えたことになる。図 A.4 の k − α σ が図A.4 に示す確率分布を持つ場合、どのあたりに目標値(どれぐらいの確率で臨界になるかというしきい値の線)を置けば良いだろうか。仮に k − α σ の確率分布の中で、左側、つまり

k − α σ を小さい側に設定したとする。これは、 k (実効増倍率の平均値)を低い側に見積もったこと、それだけ小さい実効増倍率でしか臨界にならないと評価したことに相当するため、非安全側の設定になると言える。従って、 k − α σ の分布において、右側、つまり、

k − α σ を大きい側に設定する必要がある。では、どこまで大きく設定すれば良いのだろう41

か。これは、ある確率 γ%を考えて、 k − α σ が設定した値より大きくなる確率が γ%以下になるように設定してやれば良い。このように設定したしきい値、あるいは目標値が臨界になる確率が p%の線に来るように α を調節してやれば良い。
以上の議論を踏まえると、 k sub を決定する際には以下のような手順を踏むこととなる。
①n 個の標本値より平均値 k を計算する②臨界になる確率 p を設定し、実効増倍率 k eff が平均 µ 、標準偏差 σ の正規分布に従うとして、正規分布の下側累積確率 p に対応するパーセント点 − α p を求める。
③ k − α σ ≤ µ − α p σ を満足するように、 α の値を設定することを考える。ここで、この条件式を以下のように式変形する。

k −µ≤ (α − α p ) n
σ n

(A.5)

④(A.5)式の左辺の量は、A.1 節で述べた k の統計的性質より、平均が 0、標準偏差が 1 の標

k −µの信頼区間として下側累積確率が (1 − γ )σ n準正規分布に従うことがわかる。よって、となる区間(上側累積確率が γ となる区間) α γ を、標準正規分布の上側累積確率 γ に対応するパーセント点から求める。

k −µ
≤ αγ
σ n

(A.6)

⑤(A.5), (A.6)式の右辺の値を比較することで、 k − α σ ≤ µ − α p σ を満足するためには、パラメータ α の値を以下のように設定すればよいことが分かる。

α =αp +

αγ

(A.7)

n

⑥こうして求められた α を使用して推定臨界下限増倍率を k sub (k , σ ) = k − α σ と設定する。

(A.7)式からわかるように、実効増倍率の平均値として母集団の平均 µ のかわりに標本平42

均 k を用いた場合には、 k の不確かさのぶんとして α γn だけ α の値を α p より大きく設定すれば良いことが分かる。
上 述 し た 手 順 に 基 づ い て 例 題 の 推 定 臨 界 下 限 増 倍 率 k sub (k ,σ ) を 算 出 し て み る 。

p = 2.5% 、 (1 − γ ) = 97.5% の場合には、正規分布の下/上側累積確率のパーセント点より1.96 

α p = α γ = 1.96 となるため、 k sub (k , σ ) = 1.0143 − 1.96 + × 0.02 = 0.9627 と求めら10 
れる。この値は、A.2 節で求められた k sub ( µ , σ ) の値 0.9708 より小さく、 k の不確かさも含
めて保守的に推定臨界下限増倍率を設定していると言える。手順④で(A.6)式の信頼区間を考えたことから理解できるように、もし仮に 10 個で 1 セットの臨界実験の解析を別途 1000
ケース実施したとすると、期待値として 1000 ケースのうち 25 のケース(=1000×γ＝1000×2.5%)について「設定した k sub を用いた場合、臨界になる確率が p = 2.5% を超える可能性がある」といえる。参考までに、10 個で 1 セットの臨界実験の解析を別途 1000 ケース
実施したことを想定してランダムサンプリングを行い、各ケースで k sub (k , σ ) を算出した結果を図 A.5 に示す。図 A.5 に示したように、全 1000 ケースのうち 26 ケースが A.2 で求めた「正しい」推定臨界下限増倍率である k sub ( µ , σ ) = 0.9708 よりも大きく設定される結果となり、その頻度は γ = 2.5% にほぼ等しいことが確認できる。つまり、推定臨界下限増倍率の「設定に失敗する」確率が予測した 2.5%とほぼ一致していることが分かる。

推定臨界下限増倍率

0.99
0.98
0.97
0.96
0.95
0.94
0.93
0

200

400

600

800

1000

ケ ース
図 A.5 k sub (k , σ ) の散布図(標本数 n = 10 , p = γ = 2.5% )43A.5 実効増倍率の母集団の平均値も標準偏差も不明な場合
最後に、母集団の平均値も標準偏差も不明で、k と s の両方に不確かさがある場合について考える(図 A.6)。

k −α s
確率分布

keff
頻度分布

ασ
この割合
がp%

この割合
がγ%

図 A.6 実効増倍率の母集団の平均値も標準偏差も不明な場合
2

この場合、A.4 節で述べた内容に加えて、n 個の標本から推定された不偏分散 s がカイ二乗分布に従うことも利用して、推定臨界下限増倍率 k sub (k , s ) = k − α s のパラメータ α を設定することとなる。文献(奥野, 内藤, "臨界安全安全解析コードシステム JACS の計算誤差
評価", JAERI-M 87-057)を参考として、具体的な k sub (k , s ) の設定手順を以下にまとめる。
2①n 個の標本値より k と s を求める。
②臨界になる確率 p を設定し、実効増倍率 k eff が平均 µ 、標準偏差 σ の正規分布に従うとして、正規分布の下側累積確率 p に対応するパーセント点 − α p を求める。
③ k − α s ≤ µ − α p σ を満足するように、α の値を設定することを考える。ここで、この条件式を以下のように式変形する。

k −µ

σ

n

+αp n

(n − 1) s 2
σ2

1
n −1

≤α n

(A.8)

44

④(A.8)式の左辺の分子 X =

がα p

k −µ

σ

n

+ α p n は、A.1 節で述べた k の統計的性質より、平均n 、標準偏差が 1 の正規分布に従う。また、(A.8)式の左辺の分母 Y =(n − 1) s 2 は
σ2

A.1 節で述べた s の統計的性質より、自由度 (n − 1) のカイ二乗分布に従う。従って、A.1
2

節の(3)で述べた性質(A.4)式より、(A.8)式左辺の X
ラメータが α p

Y (n − 1) は自由度が (n − 1) 、非心パ

n の非心 t 分布に従うことが分かる。よって、 X

Y (n − 1) の信頼区間として下側累積確率が (1 − γ ) となる区間 α γ を、非心 t 分布の上側累積確率 γ に対応する
パーセント点から求める。

X≤ αγ
Y (n −1)(A.9)⑤(A.8), (A.9)式の右辺の値を比較することで、 k − α s ≤ µ − α p σ を満足するためにはパラ
メータ α の値を以下のように設定すれば良いことが分かる。

α=αγ(A.10)n⑥こうして求められた α を使用して、推定臨界下限増倍率を k sub ( k , s ) = k − α s と設定する。

上述した手順に基づいて例題の推定臨界下限増倍率 k sub (k , s ) を以下で算出してみる。

p = 2.5% の場合、 α p = 1.96 となる。また標本数 n = 10 の場合、 X

Y 9 は自由度が 9、非心パラメータが 6.198 の非心 t 分布に従うので、下側累積確率 (1 − γ ) = 97.5% に対応するパーセント点が α γ = 12.0195 と求められる。よって、k sub (k , s ) の推定で用いるパラメー
タ α は α = 12.0195

10 = 3.801 と 得 ら れ る 。 最 終 的 に 、 推 定 臨 界 下 限 増 倍 率 は

k sub (k , s ) = 1.0143 − 3.801× 0.0209 = 0.9349 と求めることができる。こうして評価された
値は、A.4 節で求めた k sub (k , σ ) の値 0.9627 よりさらに小さな値であり、 k だけでなく s の不確かさを見込んだぶんさらに保守的に k sub を設定していると言える。参考までに、10 個
で 1 セットの臨界実験の解析を別途 1000 ケース実施したことを想定してランダムサンプリ
ングを行い、各ケースで k sub (k , s ) を算出した結果を図 A.7 に示す。図 A.7 に示したように、
全 1000 ケースのうち 24 ケースが A.2 で求めた「正しい」推定臨界下限増倍率である

k sub ( µ , σ ) = 0.9708 よりも大きく設定される結果となり、その頻度は γ = 2.5% にほぼ等し
45

いことが確認できる。つまり、推定臨界下限増倍率の「設定に失敗する」確率が予測した2.5%とほぼ一致していることが分かる。

推定臨界下限増倍率

1.00
0.98
0.96
0.94
0.92
0.90
0.88
0.86
0

200

400

600

800

1000

ケ ース
図 A.7 k sub (k , s ) の散布図(標本数 n = 10 , p = γ = 2.5% )46