アルゴリズム(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の計算手法で、Unified equations for the slope, intercept, and standard error of the best straight lineに述べられています。
  • FV法
    Fitting straight lines with errors on both coordinates で述べているG. Fasano and R. Vio の計算方法を使います。
  • Deming法
    Deming回帰法は変数誤差の最尤推定法で、X/Y 誤差は、独立で同一分布にしたがっていると仮定します。
  • XおよびYエラーの相関
    XおよびYエラーの相関r_i (York法のみ)
  • X/Yの標準偏差
    X/Yの標準偏差(Deming 法のみ)

値(York 法)

線形フィットを実行すると、分析レポートシートに計算された値が出力されます。 パラメータ表には、モデルの傾きと切片(括弧内の数字は生成された値を示す)が表示されます。

フィットパラメータ

York Error.png

フィット値と標準誤差

xとyの重み(誤差)に関するW_iを定義します。

W_i = \frac{\sigma x_i\sigma y_i}{\sigma x_i+\beta_1^2\sigma y_i-2\beta_1 r_ia_i}

(3)

このとき、r_iは、XとYの誤差の相関係数 (すなわち、\sigma x_i\sigma y_i)、および、a_i=\sqrt{\sigma x_i\sigma y_i}です。

重み付け(誤差)無しの(X_i, Y_i)のフィット線の傾きは、\beta_1\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_ib_iV_i}}{\sum{W_ib_iU_i}}

(5)

U および V はX および Yの偏差です。

\left\{\begin{matrix}
U=X-\hat{X}\\
V=Y-\hat{Y}
\end{matrix}\right.

および、

b_i=W_i[\frac{U_i}{\sigma_{y_i}}+\frac{\hat{\beta_1}}{\sigma_{x_i}}{V_i}-(\beta U_i+V_i)\frac{r_i}{a_i}]

(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_iU_i^2}}

(8)

パラメータとしての標準誤差は、次のように与えられます。

\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-bX_i-a)^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\,\!

The t-values can be computed by:

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)

ここでUCLLCLは、それぞれ上側信頼区間下側信頼区間のことです。

CI 半幅

信頼区間の半値幅は以下の通りです。

CI=\frac{UCL-LCL}2

(16)

ここでUCLとLCLは、それぞれ上側信頼区間下側信頼区間です。

詳細については、以下の参考文献 1 をご覧ください。

フィット統計

York Stats.png

自由度

df=n-2

(17)

nは、ポイントの全数です。

残差平方和

RSS=\sum^n_{i=1} \frac{(\beta_0+\beta_1x_i-y_i)^2}{\sigma^2_{yi}+\beta_1^2\sigma^2_{xi}}

(18)

自由度あたりカイ二乗

\sigma^2=\frac{RSS}{n-2}

(19)

ピアソンのr

単純な線形回帰では、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法)

Fitting straight lines with errors on both coordinates で述べているG. Fasano and R. Vio の計算方法を使います。

重みは、次のように定義されます。

W_i=\frac{1}{\beta_1^2\sigma {x_{i}}^2+\sigma {x_{i}}^2}

(25)

重み付け(誤差)無しの(X_i, Y_i)のフィット線の傾きは、\beta_1です。

次のようにします。

\bar{x}=\frac{\sum{W_ix_i}}{\sum W_i}

(26)

\bar{y}=\frac{\sum{W_iy_i}}{\sum W_i}

(27)

合計K^2=\sum{W_i (y_i-\beta_0-\beta_1 x_i)}を最小化し、部分微分を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 をご覧ください。

値 (Deming法)

線形フィットを実行すると、分析レポートシートに計算された値が出力されます。 パラメータ表には、モデルの傾きと切片(括弧内の数字は生成された値を示す)が表示されます。

フィットパラメータ

Deming Error.png

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~ N(0,\sigma^2)と独立で同一の分布に従い、\sigma y_iは、(0,\lambda \sigma^2)独立で同一の分布に従います。\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 and 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)

ここでUCLLCLは、それぞれ上側信頼区間下側信頼区間のことです。

CI 半幅

信頼区間の半値幅は以下の通りです。

CI=\frac{UCL-LCL}2

(41)

ここでUCLとLCLは、それぞれ上側信頼区間下側信頼区間です。

詳細については、以下の参考文献 1 をご覧ください。

フィット統計

Deming Stats.png

自由度

df=n-2

(42)

nは、ポイントの全数です。

残差平方和

式(33)を参照。

自由度あたりカイ二乗

\sigma^2=\frac{RSS}{n-2}

(43)

ピアソンのr

単純な線形回帰では、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

Covariance and Correlation Matrix

線形回帰の共分散行列は次のように計算されます。


\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検索

残差プロット

残差と独立変数

残差散布図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 はデータセットの合計数で、ii 番目のデータです。 なお、正規確率プロットとQ-Qプロットについてをご覧ください。

参考文献

  1. 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).
  2. 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.