数値計算の精度

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: Statistics タブに表示される mean と std を比較する
   • Linear Regression: Linear Regression タブで Response Variable（応答変数）に y を、Predictor Variables（説明変数）に x を設定し、係数や R-squared を比較する。Filip (Orthogonal Poly.) 行については、Orthogonal Polynomials タブで x から 10 次の直交多項式列を生成し、それを Predictor Variables として使用する。R² と RMSE は基底によらず不変なので UI 上の値を直接比較できます（F 統計量は UI に表示されないため比較できません）。係数と SE は直交基底での値が表示されるため、NIST 認定値（生多項式基底）とは直接比較できません。表の LRE(係数) と LRE(SE) は自動テストで算出した値です
   • ANOVA: ANOVA タブで group 列を Factor A（因子）に、value 列を Response Variable（応答変数）に設定し、ANOVA Table に表示される SS を比較する

精度の読み方

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

$$\text{LRE} = -\log_{10} \frac{|\text{計算値} - \text{認定値}|}{|\text{認定値}|}$$

認定値が 0 の場合や統計量が定義できない場合はこの式を適用できません。表中の "exact" は、NIST 認定値が 0 であるか統計量が定義できないため、LRE を計算できないことを示します。

MIDAS を含むすべてのブラウザアプリケーションは IEEE 754 倍精度浮動小数点数で計算します。仮数部は 52 ビットを格納し、暗黙の先頭ビットを含めた有効ビット数は 53 です。LRE の理論上の上限は約 15.9 です。表の LRE は小数第 1 位に丸めているため、15.95 以上の値は 16.0 と表示されます。なお相対誤差が厳密に 0（完全一致）の場合、上式は発散するため、MIDAS の算出ルーチンは便宜上の有限値として LRE = 15 を返します。

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

基本統計量

MIDAS の Statistics タブに表示される mean と std を NIST 認定値と比較しました。標準偏差は標本標準偏差で $n - 1$ で割ります。

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

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

線形回帰データセット

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

Wampler1 と Wampler2 はノイズなしのデータで、モデルがデータに完全に一致します。表中の "exact" は 2 つの状況を表します。LRE(SE) と LRE(RMSE) の "exact" は、SE と RMSE の NIST 認定値が 0 のため LRE を計算できないことを示します。LRE(F) の "exact" は、MSE = 0 のため F 統計量が定義できないことを示します。

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

一元配置分散分析

MIDAS の ANOVA タブで各データセットの一元配置分散分析を実行し、得られた群間平方和（SS_B）と群内平方和（SS_W）を NIST 認定値と比較しました。表の数値は LRE です。「難易度」は NIST StRD による計算難易度の分類で、データ値のオフセットの大きさに対応します。

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

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

既知の制約

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

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

SmLs07-SmLs09 (ANOVA Higher): データ値が $10^{12}$ 台で、変動は小数点以下のみです。群間平方和（SS_B）は間接法（SS_total − SS_within）で計算します。SmLs07-09 では SS_total と SS_within がほぼ等しいため、その差分計算で桁落ちが発生します。SS_total と SS_within の各算出には Welford のオンラインアルゴリズム を使用しており、大オフセットの影響を軽減していますが、差分での桁落ちは避けられません。SS_B の精度は約 4 有効桁です。群内平方和（SS_W）は各群内の Welford 分散から積算するため、大オフセットの影響を軽減して約 4 有効桁が得られます。

データの出典

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

このページの Markdown 版もあります。