アルゴリズム: Xエラーあり線形フィット
フィッティングモデル
与えられたデータセット\((X_i,Y_i), (\sigma_{x_i},\sigma_{y_i}), i=1,2,\ldots n\)において、Xは独立変数、Yは従属変数であり、\((\sigma_{x_i},\sigma_{y_i})\)はX,Yそれぞれのエラーを表します。 Xエラー付き線形フィットは、データを次の形式のモデルに適合させます。
| \[y=\beta _0+\beta _1x+\varepsilon\] |
(1) |
|---|---|
| \[\left\{\begin{matrix} x_i=X_i+\sigma_{x_i}\\ y_i=Y_i+\sigma_{y_i} \end{matrix}\right.\] |
(2) |
フィット制御
計算法
- York法
- York法はD.Yorkの計算法であり、最良の傾き、切片、標準偏差を持つ最善の直線について説明されています。
- FV法
- FV法はGiocanniFasano & Roberto Vioによる計算方法で、両方の座標にエラーのある直線へのフィットについて説明されています。
- Deming法
- Deming回帰は、変数内誤差モデルの最尤推定であり、X/Y 誤差は独立して等分布していると仮定されます。
- XおよびYエラーの相関
- XおよびYエラーの相関\(r_i\)(York法のみ)
- X/Yの標準偏差
- X/Yの標準偏差(Deming法のみ）
値(York法)
線形フィットを実行すると、分析レポートシートに計算された値が出力されます。 パラメータ表にはモデルの傾きと切片（カッコ内の数字は計算された値を示す）が表示されます。
フィットパラメータ
フィット値と標準誤差
xとyの重み（誤差）に関する\(W_i\)を定義します。
| \[W_i = \frac{\omega_{x_i}\omega_{y_i}}{\omega_{x_i}+\beta_1^2\omega_{y_i}-2\beta_1 r_ia_i} =\frac{1}{\sigma_{y_i}^2+\beta_1^2\sigma_{x_i}^2 - 2\beta_1 r_i \sigma_{x_i} \sigma_{y_i}}\] |
(3) |
|---|
このとき\(\omega_{x_i}=\frac{1}{\sigma_{x_i}^2}, \ \omega_{y_i}=\frac{1}{\sigma_{y_i}^2}\)は\((X_i, Y_i)\)の重みで、\(r_i\)はXとYの誤差の相関係数(すなわち\(\sigma_{x_i}\)と\(\sigma_{y_i}\))で\(\alpha_i=\sqrt{\omega_{x_i} \omega_{y_i}}\)です。
重みづけ（誤差）のない\((X_i, Y_i)\)のフィット線の傾きは\(\beta_1\)の初期値になります。設定した許容値\(\beta_1\)に収まるまでこれらは反復計算されます。
X_Y誤差のあるもっとフィットする線の推定パラメータ\(\hat{\beta_0}\)と\(\hat{\beta_1}\)の簡潔な方程式は次のようになります。
| \[\hat{\beta_0}=\bar{Y}-\hat{\beta_1}\bar{X}\] |
(4) |
|---|
| \[\hat{\beta_1}=\frac{\sum{W_i b_i V_i}}{\sum{W_i b_i U_i}}\] |
(5) |
|---|
ここで\(\bar{X} = \frac{ \sum{W_i X_i} }{ \sum{W_i} }, \ \bar{Y} = \frac{ \sum{W_i Y_i} }{ \sum{Y_i} }\)です。
UとVはXとYの偏差です。
\[\left\{\begin{matrix} U=X-\bar{X}\\ V=Y-\bar{Y} \end{matrix}\right. \]
および、
| \[b_i=W_i \left[\frac{U_i}{\omega_{y_i}}+\frac{\hat{\beta_1}}{\omega_{x_i}}{V_i}-(\beta U_i+V_i)\frac{r_i}{\alpha_i} \right]\] |
(6) |
|---|
パラメータの対応する変数\(\sigma^2\)と標準誤差\(s\)は次のようになります。
| \[\sigma_{\hat{\beta_0}}^2=\frac{1}{\sum{W_i}}+\bar{x}^2\sigma_{\hat{\beta_1}}^2\] |
(7) |
|---|---|
| \[\sigma_{\hat{\beta_1}}^2=\frac{1}{\sum{W_i u_i^2}}\] |
(8) |
ここで\(\bar{x} = \frac{ \sum{W_i x_i} }{ \sum{W_i} }\)、\(x_i\)は\(X_i\)の推定値で、\(u_i=x_i - \bar{x}\)です。
パラメータの標準誤差は最終的に次のようになります。
| \[\varepsilon_{\hat{\beta_0}}=\sqrt{\sigma _{\hat{\beta_0}}^2}\sqrt{\frac{S}{n-2}}\] |
(9) |
|---|
| \[\varepsilon_{\hat{\beta_1}}=\sqrt{\sigma _{\hat{\beta_1}}^2}\sqrt{\frac{S}{n-2}}\] |
(10) |
|---|
ここで\(S\)は
| \[S=\sum W_i(Y_i - \beta_1 X_i- \beta_0)^2\] |
(11) |
|---|
t値と信頼水準
回帰の仮定が成り立つ場合次のようになります。
| \(\frac{{\hat \beta _0}-\beta _0}{\varepsilon _{\hat \beta _0}}\sim t_{n^{*}-1}\) かつ \(\frac{{\hat \beta _1}-\beta _1}{\varepsilon _{\hat \beta _1}}\sim t_{n^{*}-1}\) |
(12) |
|---|
フィッティングパラメータが0ではないことを調べるためにt 検定を使うことができます。これは、 \(\beta _0= 0\,\!\) (真ならば、フィット直線が原点を通る) または \(\beta _1= 0\,\!\) であるかどうかを検定します。t 検定の仮説検定は次のようになります。
- \(H_0 : \beta _0= 0\,\! \) \(H_0 : \beta _1= 0\,\!\)
- \(H_\alpha : \beta _0 \neq 0\,\!\) \(H_\alpha : \beta _1 \neq 0\,\!\)
t 値は、次の式で計算できます。
| \(t_{\hat \beta _0}=\frac{{\hat \beta _0}-0}{\varepsilon _{\hat \beta _0}}\) かつ \(t_{\hat \beta _1}=\frac{{\hat \beta _1}-0}{\varepsilon _{\hat \beta _1}}\) |
(13) |
|---|
計算されたt値を使って対応する帰無仮説を棄却するかどうかを決定できます。通常与えられた有意水準\(\alpha\,\!\)に対して、\(|t|>t_{\frac \alpha 2}\)のときに\(H_0 \,\!\)を棄却します。また、p値または有意水準がt検定とともに出力されます。p値が\(\alpha\,\!\)よりも小さい場合に帰無仮説\(H_0 \,\!\)を棄却することができます。
Prob>|t|
上記のt検定における\(H_0 \,\!\)が真である確率
| \[prob=2(1-tcdf(|t|,df_{Error}))\,\!\] |
(14) |
|---|
ここで、tcdf(t, df) は、自由度 df を持つスチューデントt 分布の下側の確率を計算します。
LCLとUCL
t値から各パラメータの信頼区間\((1-\alpha )\times 100\%\)を次の式で計算できます。
| \[\hat \beta _j-t_{(\frac \alpha 2,n^{*}-k)}\varepsilon _{\hat \beta _j}\leq \hat \beta _j\leq \hat \beta _j+t_{(\frac \alpha 2,n^{*}-k)}\varepsilon _{\hat \beta _j}\] |
(15) |
|---|
ここで\(UCL\)と\(LCL\)はそれぞれ上限信頼区間と下限信頼区間の略です。
CI 半幅
信頼区間の半値幅は以下の通りです。
| \[CI=\frac{UCL-LCL}2\] |
(16) |
|---|
ここでUCLとLCLは、それぞれ上側信頼区間と下側信頼区間です。
より詳細は参考文献1（下記）をご覧ください。
フィット統計
自由度
| \[df=n-2\] |
(17) |
|---|
nは合計ポイント数です。
残差平方和
| \[RSS=\sum^n_{i=1} \frac{(\beta_0+\beta_1 x_i - y_i)^2}{\sigma^2_{y_i}+\beta_1^2\sigma^2_{x_i}}\] |
(18) |
|---|
自由度あたりカイ二乗
| \[\sigma^2=\frac{RSS}{n-2}\] |
(19) |
|---|
ピアソンのｒ
単純な線形回帰では、xとyの相関係数は、r で表され、次の式に等しくなります。
| \(r=R\,\!\) \(\beta _1\,\!\) が正の場合 |
(20) |
|---|---|
| \(r=-R\,\!\) \(\beta _1\,\!\) が負の場合 |
\[R^2\] これは次式のように計算されます。
| \[R^2=\frac{SXY}{SXX*TSS}=1-\frac{RSS}{TSS}\] |
(21) |
|---|---|
| \[TSS=\sum_{i=1}^n(y_i-\bar{y})^2\] |
Root MSE(SD)
誤差の平均平方の平方根または、残差標準偏差は、次式に等しくなります。
| \[RootMSE=\sqrt{\frac{RSS}{df_{Error}}}\] |
(22) |
|---|---|
| \[df_{Error}=n-2\] |
共分散行列と相関行列
線形回帰における共分散行列は次のように計算されます。
| \[ \begin{pmatrix} Cov(\beta _0,\beta _0) & Cov(\beta _0,\beta _1)\\ Cov(\beta _1,\beta _0) & Cov(\beta _1,\beta _1) \end{pmatrix}=\sigma ^2\frac 1{SXX}\begin{pmatrix} \sum \frac{x_i^2}n & -\bar x \\-\bar x & 1 \end{pmatrix}\] |
(23) |
|---|
2つのパラメータ間の相関は、
| \[ \rho (\beta _i,\beta _j)=\frac{Cov(\beta _i,\beta _j)}{\sqrt{Cov(\beta _i,\beta _i)}\sqrt{Cov(\beta _j,\beta _j)}} \] |
(24) |
|---|
値 (FV法)
FV法はGiocanniFasano & Roberto Vioによる計算方法で、両方の座標にエラーのある直線へのフィットについて説明されています。
重みは次のように定義されます。
| \[W_i=\frac{1}{\beta_1^2\sigma_{x_{i}}^2+\sigma_{y_{i}}^2}\] |
(25) |
|---|
重みづけ（誤差）のない\((X_i, Y_i)\)のフィット線の傾きは\(\beta_1\)です。
次のようにします。
| \[\bar{x}=\frac{\sum{W_i x_i}}{\sum W_i}\] |
(26) |
|---|---|
| \[\bar{y}=\frac{\sum{W_i y_i}}{\sum W_i}\] |
(27) |
合計\(K^2=\sum{W_i (y_i-\beta_0-\beta_1 x_i)^2}\)を最小化し、偏微分を0に設定することで値\(\beta_0\)と\(\beta_1\)を推定することができます。
| \[\hat{\beta_0}=\bar{y}-\hat{\beta_1}\bar{x}\] |
(28) |
|---|---|
| \[a\hat{\beta_1}^2+b\hat{\beta_1}-c=0\] |
(29) |
ここで
| \[a=\sum{W_i^2\sigma_{x_i}^2(y_i-\bar{y_i})(x_i-\bar{x_i})}\] |
(30) |
|---|---|
| \[b=\sum{W_i^2[\sigma_{y_i}^2(x_i-\bar{x_i})^2-\sigma_{x_i}^2(y_i-\bar{y_i})^2]}\] |
(31) |
| \[c=\sum{W_i^2\sigma_{y_i}^2(y_i-\bar{y_i})(x_i-\bar{x_i})}\] |
(32) |
設定した許容値\(\hat{\beta_1}\)に収まるまで\(\hat{\beta_1}\)は反復計算されます。
それぞれのパラメータの標準誤差については線形回帰モデルを参照してください。
より詳細は参考文献2（下記）をご覧ください。
5値 (Deming法)
線形フィットを実行すると、分析レポートシートに計算された値が出力されます。 パラメータ表にはモデルの傾きと切片（カッコ内の数字は計算された値を示す）が表示されます。
フィットパラメータ
Deming回帰はxとyに測定誤差あることが前提とされる場合に使われます。
| \[y=\beta _0+\beta _1x+\varepsilon\] | |
|---|---|
| \[\left\{\begin{matrix} x_i=X_i+\sigma_{x_i}\\ y_i=Y_i+\sigma_{y_i} \end{matrix}\right.\] |
\(\sigma_{x_i}\)は、\(\sigma_{x_i} \sim \mathcal{N}(0,\sigma^2)\)と独立で同一の分布に従い、\(\sigma_{y_i}\)は、\(\sigma_{y_i} \sim \mathcal{N}(0,\lambda \sigma^2)\)と独立で同一の分布に従うと仮定します。 ここで、\(\mathcal{N}(0,\sigma^2)\)は平均0と標準偏差\(\sigma\)の正規分布を示します。\(\lambda=1\)の場合それは直交回帰です。 モデルの加重残差平方和は最小化され、
| \[RSS=\sum^n_{i=1}\left ((x_i-X_i)^2+\frac{(y_i-\beta_0-\beta_1X_i)^2}{\lambda}\right)\] |
(33) |
|---|
フィット値と標準誤差
パラメータを解くことができます。
| \[\hat{\beta_1}=\frac{SYY-\lambda SXX+\sqrt{(SYY-\lambda SXX)^2+4\lambda SXY^2}}{2SXY}\] |
(34) |
|---|
| \[\hat{\beta_0}=\bar{y}-\hat{\beta_1}\bar{x}\] |
(35) |
|---|
ここで、
| \[\bar{x}=\frac{1}{n}\sum_{i=1}^2{x_i}, \bar{y}=\frac{1}{n}\sum_{i=1}^n{y_i}\] |
|---|
| \[u_i=x_i-\bar{x}\] | |
|---|---|
| \[v_i=y_i-\bar{y}\] |
さらに
| \[SXX=\sum_{i=1}^n u_i^2\] | |
|---|---|
| \[SYY=\sum_{i=1}^n v_i^2\] | |
| \[SXY=\sum_{i=1}^n u_iv_i\] |
パラメータの対応する変数は次のようになります。
| \[\sigma^2_{\hat \beta _0}=\frac{1}{nw}+2(\bar{x}+2\bar{z})\bar{z}Q+(\bar{x}+2\bar{z})^2 \sigma_{\bar{\beta_1}}^2\] | |
|---|---|
| \[\sigma^2_{\hat \beta _1}=Q^2w^2\sigma^2\sum^n_{i=1}(\lambda u_i^2+v_i^2)\] |
パラメータの標準誤差は次のように推定できます。
| \[\varepsilon _{\hat \beta _0}=\sqrt{\sigma^2_{\hat \beta _0}}\] |
(37) |
|---|---|
| \[\varepsilon _{\hat \beta _1}=\sqrt{\sigma^2_{\hat \beta _1}}\] |
(38) |
および、
| \[w=\frac{1}{\sigma^2(\lambda+\hat{\beta_1}^2)}\] |
|---|
| \[z_i=w\sigma^2(\lambda u_i+\hat{\beta_1} v_i)\] |
|---|
| \[\bar{z}=\frac{1}{n}\sum_{i=1}^n z_i\] |
|---|
| \[Q=\frac{1}{w \sum_{i=1}^n \left(\frac{u_iv_i}{\hat{\beta_1}}+4(z_i-\bar{z})(z_i-u_i)\right)}\] |
|---|
| \[\sigma=\sqrt{\frac{\sum^n_{i=1}(x_i-X_i)^2+\frac{\sum^n_{i=1}(y_i-\hat{\beta_0}-\hat{\beta_1}x_i)^2}{\lambda}}{n-2}}\] |
|---|
t値と信頼水準
回帰の仮定が成り立つ場合次のようになります。
| \(\frac{{\hat \beta _0}-\beta _0}{\varepsilon _{\hat \beta _0}}\sim t_{n^{*}-1}\) かつ \(\frac{{\hat \beta _1}-\beta _1}{\varepsilon _{\hat \beta _1}}\sim t_{n^{*}-1}\) |
フィッティングパラメータが0ではないことを調べるためにt 検定を使うことができます。これは、 \(\beta _0= 0\,\!\) (真ならば、フィット直線が原点を通る) または \(\beta _1= 0\,\!\) であるかどうかを検定します。t 検定の仮説検定は次のようになります。
- \(H_0 : \beta _0= 0\,\! \) \(H_0 : \beta _1= 0\,\!\)
- \(H_\alpha : \beta _0 \neq 0\,\!\) \(H_\alpha : \beta _1 \neq 0\,\!\)
t 値は、次の式で計算できます。
| \(t_{\hat \beta _0}=\frac{{\hat \beta _0}-0}{\varepsilon _{\hat \beta _0}}\) かつ \(t_{\hat \beta _1}=\frac{{\hat \beta _1}-0}{\varepsilon _{\hat \beta _1}}\) |
(38) |
|---|
計算されたt値を使って対応する帰無仮説を棄却するかどうかを決定できます。通常与えられた有意水準\(\alpha\,\!\)に対して、\(|t|>t_{\frac \alpha 2}\)のときに\(H_0 \,\!\)を棄却します。また、p値または有意水準がt検定とともに出力されます。p値が\(\alpha\,\!\)よりも小さい場合に帰無仮説\(H_0 \,\!\)を棄却することができます。
Prob>|t|
上記のt検定における\(H_0 \,\!\)が真である確率
| \[prob=2(1-tcdf(|t|,df_{Error}))\,\!\] |
(39) |
|---|
ここで、tcdf(t, df) は、自由度 df を持つスチューデントt 分布の下側の確率を計算します。
LCLとUCL
t値から各パラメータの信頼区間\((1-\alpha )\times 100\%\)を次の式で計算できます。
| \[\hat \beta _j-t_{(\frac \alpha 2,n^{*}-k)}\varepsilon _{\hat \beta _j}\leq \hat \beta _j\leq \hat \beta _j+t_{(\frac \alpha 2,n^{*}-k)}\varepsilon _{\hat \beta _j}\] |
(40) |
|---|
ここで\(UCL\)と\(LCL\)はそれぞれ上限信頼区間と下限信頼区間の略です。
CI 半幅
信頼区間の半値幅は以下の通りです。
| \[CI=\frac{UCL-LCL}2\] |
(41) |
|---|
ここでUCLとLCLは、それぞれ上側信頼区間と下側信頼区間です。
より詳細は参考文献1（下記）をご覧ください。
フィット統計
自由度
| \[df=n-2\] |
(42) |
|---|
nは合計ポイント数です。
残差平方和
式(33)を参照
自由度あたりカイ二乗
| \[\sigma^2=\frac{RSS}{n-2}\] |
(43) |
|---|
ピアソンのｒ
単純な線形回帰では、xとyの相関係数は、r で表され、次の式に等しくなります。
| \(r=R\,\!\) \(\beta _1\,\!\) が正の場合 |
(44) |
|---|---|
| \(r=-R\,\!\) \(\beta _1\,\!\) が負の場合 |
\[R^2\] これは次式のように計算されます。
| \[R^2=\frac{SXY}{SXX*TSS}=1-\frac{RSS}{TSS}\] |
(45) |
|---|---|
| \[TSS=\sum_{i=1}^n(y_i-\bar{y})^2\] |
Root MSE(SD)
誤差の平均平方の平方根は、次式に等しくなります。
| \[RootMSE=\sqrt{\frac{RSS}{df_{Error}}}\] |
(46) |
|---|---|
| \[df_{Error}=n-2\] |
共分散行列と相関行列
線形回帰における共分散行列は次のように計算されます。
| \[ \begin{pmatrix} Cov(\beta _0,\beta _0) & Cov(\beta _0,\beta _1)\\ Cov(\beta _1,\beta _0) & Cov(\beta _1,\beta _1) \end{pmatrix}=\begin{pmatrix} \ \sigma^2_{\hat{\beta_0}} & -\bar{x}\sigma^2_{\hat \beta _1} \\-\bar{x}\sigma^2_{\hat \beta _1} &\sigma^2_{\hat{\beta_1}} \end{pmatrix}\] |
(47) |
|---|
2つのパラメータ間の相関は、
| \[ \rho (\beta _i,\beta _j)=\frac{Cov(\beta _i,\beta _j)}{\sqrt{Cov(\beta _i,\beta _i)}\sqrt{Cov(\beta _j,\beta _j)}} \] |
(48) |
|---|
X/Y検索
残差プロット
残差vs.独立
残差\(res\) vs. 独立変数\(x_1,x_2,\dots,x_k\)の散布図プロットではそれぞれのプロットは別々のグラフに配置されます。
残差vs.予測値
残差\(res\)の散布図プロット vs. フィット結果\(\hat{y_i}\)
残差vs.データ順序
\(res_i\) vs. 順序\(i\)
残差のヒストグラム
残差\(res_i\)のヒストグラムプロット
残差のラグプロット
残差\(res_i\) vs. ラグ残差\(res_{(i–1)}\)
正規残差確率プロット
残差の正規確率プロットは分散が成否分布しているかどうかを調べるために使います。結果のプロットはおおよそ線形で、誤差範囲は正規分布していると仮定することができます。プロットはパーセンタイル対順序化された残差をベースにしており、パーセンタイルは次のように仮定されます。
\[\frac{(i-\frac{3}{8})}{(n+\frac{1}{4})}\]
ここで、n はデータセットの合計数で、iはi番目のデータを表します。 確率プロットとQ-Qプロットもご覧ください。
参考文献
- York D, Unified equations for the slope, intercept, and standard error of the best straight line, American Journal of Physics, Volume 72, Issue 3, pp. 367-375 (2004).
- G. Fasano and R. Vio, "Fitting straight lines with errors on both coordinates", Newsletter of Working Group for Modern Astronomical Methodology, No. 7, 2-7, Sept. 1988.