Language

数値計算の精度

MIDAS の統計計算が正確かどうかを、ユーザー自身で確認する方法を説明します。

NIST Statistical Reference Datasets

NIST Statistical Reference Datasets (StRD) は米国国立標準技術研究所が統計ソフトウェアの精度検証用に公開しているベンチマークデータセットです。各データセットに 15 有効桁の認定値が付属しています。

NIST StRD には Univariate、Linear Regression、Nonlinear Regression、ANOVA、MCMC の 5 カテゴリがあります。このページでは MIDAS の機能に対応する Univariate、Linear Regression、ANOVA を扱います。Nonlinear Regression と MCMC に対応する機能は MIDAS にありません。

NIST データセットに限らず、R や Python で計算した結果がわかっているデータセットでも同じ方法で検証できます。

検証方法

UI で検証する

  1. 下の表から CSV ファイルをダウンロードする
  2. MIDAS で CSV を読み込む
  3. 検証したいカテゴリに応じて操作する
    • Univariate: Data Table パネルに表示される平均と標準偏差を比較する
    • Linear Regression: Linear Regression タブで応答変数 y と説明変数 x を設定し、係数や R-squared を比較する
    • ANOVA: ANOVA タブで group 列を因子、value 列を応答変数に設定し、SS と F 統計量を比較する

精度の読み方

表の数値は Log Relative Error (LRE) です。MIDAS の計算値と NIST 認定値の相対誤差の常用対数の符号を反転した値で、一致する有効数字の桁数に相当します。

LRE=−log⁡10∣計算値−認定値∣∣認定値∣\text{LRE} = -\log_{10} \frac{|\text{計算値} - \text{認定値}|}{|\text{認定値}|}

認定値が 0 の場合や統計量が定義できない場合はこの式を適用できません。表中の "exact" は計算値と認定値がともに 0 であるか、統計量が定義できず計算値・認定値の双方が欠損であることを示します。

MIDAS を含むすべてのブラウザアプリケーションは IEEE 754 倍精度浮動小数点数で計算します。仮数部は 52 ビットを格納し、暗黙の先頭ビットを含めた有効ビット数は 53 です。LRE の理論上の上限は約 15.9 です。表の LRE は小数第 1 位に丸めているため、15.95 以上の値は 16.0 と表示されます。

この表の LRE は MIDAS の計算エンジンから算出し、NIST 認定値との比較を自動テストで継続的に検証しています。

基本統計量

MIDAS の Data Table パネルに表示される平均と標準偏差を NIST 認定値と比較しました。標準偏差は標本標準偏差で n−1n - 1 で割ります。

データセットnLRE(平均)LRE(SD)
PiDigits50001514.9
Lottery21815.215.7
Lew2001515.2
Mavro501513.1
Michelso1001513.9
NumAcc131515
NumAcc210011514.2
NumAcc3100115.99.5
NumAcc4100115.78.3

NumAcc3 と NumAcc4 は平均が 10610^6 ~ 10710^7 で標準偏差が 0.1 のデータセットです。各値から平均を引く計算で、大きさの近い数同士の減算により有効桁数が失われる桁落ちが発生し、標準偏差の精度が低下します。

各データセットの詳細は NIST StRD Univariate で公開されています。

線形回帰データセット

MIDAS で各データセットの回帰を実行し、得られた係数・標準誤差・R-squared・残差 SD・F 統計量を NIST 認定値と比較しました。表の数値は LRE です。複数の回帰係数があるデータセットでは、各係数の LRE のうち最小値を示しています。LRE(SE) も同様です。最小値を採用するのは、最も精度が低い推定値を基準に評価するためです。R²、残差 SD、F 統計量はモデル全体で 1 つの値のため、そのまま表示しています。

データセットnLRE(係数)LRE(SE)LRE(R²)LRE(残差SD)LRE(F)
Norris3612.313.815.513.911.5
Pontius4011.91316139.5
NoInt11114.715.415.715.313.9
NoInt2315.315.81615.514.2
Filip827.37.510.48.27.9
Longley161312.314.312.312
Wampler1219.5exact15exactexact
Wampler22112.6exact15exactexact
Wampler3219.513.61614.411
Wampler4217.813.515.914.815.7
Wampler5215.813.513.714.813.7

Wampler1 と Wampler2 はノイズなしのデータで、モデルがデータに完全に一致します。残差がすべて 0 のため SE と残差 SD は 0 になり、計算値と認定値が一致します。MSE = 0 のため F 統計量は定義できず、MIDAS は null を返します。表中の "exact" はこれらの完全一致・計算不能を示します。

各データセットのモデル仕様・認定値・データの説明は NIST StRD Linear Regression で公開されています。

一元配置分散分析

MIDAS の ANOVA タブで各データセットの一元配置分散分析を実行し、得られた群間平方和（SS_B）、群内平方和（SS_W）、F 統計量を NIST 認定値と比較しました。表の数値は LRE です。

データセット難易度n群数LRE(SS_B)LRE(SS_W)LRE(F)
SiRstvLower25513.413.112.9
SmLs01Lower189915.615.215
SmLs02Lower1809915.415.214.9
SmLs03Lower18009915.315.515.9
AtmWtAgAverage4821110.911.7
SmLs04Average18999.310.39.3
SmLs05Average180999.310.39.3
SmLs06Average1800999.310.39.3
SmLs07Higher18993.34.33.3
SmLs08Higher180993.34.33.3
SmLs09Higher1800993.34.33.3

SiRstv はシリコンの比抵抗の実測データ（5 群各 5 観測）です。AtmWtAg は銀の原子量測定データ（2 群、計 48 観測）です。SmLs01-SmLs09 は同一構造の生成データに異なる定数オフセットを加えたもので、定数先頭桁の数（3 桁 / 7 桁 / 13 桁）によって計算難易度が上がります。

各データセットの詳細は NIST StRD ANOVA で公開されています。

既知の制約

Filip: 10次多項式で、計画行列の条件数 κ(X′X)\kappa(X'X) が約 101410^{14}（κ(X)≈107\kappa(X) \approx 10^7）です。条件数が大きいほど丸め誤差の影響が大きくなります。素の多項式基底での係数の精度は 7 有効桁です。Orthogonal Polynomials で直交多項式列を生成し、それを説明変数に使うと条件数が約 1 に低下し、係数の精度は 10 有効桁以上に改善されます。

Wampler5: Wampler1-5 はすべて同じデザイン行列を使うため条件数は同一です。Wampler5 はノイズが最も大きく、信号対雑音比が低いため、条件数に起因する丸め誤差の影響が係数の推定値に対して相対的に大きくなります。係数の精度は 6 有効桁です。

SmLs07-SmLs09 (ANOVA Higher): データ値が 101210^{12} 台で、変動は小数点以下のみです。13 桁の定数先頭桁があるため、群平均と全体平均の差 yˉi−yˉ\bar{y}_i - \bar{y} の計算で桁落ちが発生します。群間平方和（SS_B）と F 統計量の精度は約 3 有効桁にとどまります。群内平方和（SS_W）は各群内の値と群平均の差から Welford のオンラインアルゴリズムで計算するため、大オフセットの影響を受けず約 4 有効桁が得られます。

データの出典

National Institute of Standards and Technology. (1999). Statistical Reference Datasets. Standard Reference Database 140. https://doi.org/10.18434/T43G6C