単回帰モデル ¶ ある変数x 2 x_2 x 2 を別の変数x 1 x_1 x 1 から予測する問題を考える。
x 2 = α x 1 + β + e x_2 = \alpha x_1 + \beta + e x 2 = α x 1 + β + e x 1 x_1 x 1 は予測変数 (predictor variable)、x 2 x_2 x 2 は基準変数 (criterion variable)と呼ばれる確率変数である。
e e e は予測誤差を表現する確率変数であり誤差変数 といい、以下の性質を仮定する。
パラメータの推定には最尤法や最小二乗法やモーメント法などが使われる
モーメント法による単回帰モデルの推定 ¶ 単回帰モデルの両辺の期待値をとると
E [ x 2 ] = α E [ x 1 ] + E [ β ] + E [ e ] ⏟ = 0 ⟹ β = E [ x 2 ] − α E [ x 1 ] E[x_2] = \alpha E[x_1] + E[\beta] + \underbrace{ E[e] }_{=0}\\
\implies \beta = E[x_2] - \alpha E[x_1] E [ x 2 ] = α E [ x 1 ] + E [ β ] + = 0 E [ e ] ⟹ β = E [ x 2 ] − α E [ x 1 ] となる。
次に、単回帰モデルの両辺にx 1 x_1 x 1 をかけてから期待値をとると
E [ x 2 x 1 ] = α E [ x 1 x 1 ] + β E [ x 1 ] + E [ e x 1 ] ⏟ = 0 E[x_2 x_1] = \alpha E[x_1 x_1] + \beta E[x_1] + \underbrace{ E[e x_1] }_{=0}\\ E [ x 2 x 1 ] = α E [ x 1 x 1 ] + βE [ x 1 ] + = 0 E [ e x 1 ] β = E [ x 2 ] − α E [ x 1 ] \beta = E[x_2] - \alpha E[x_1] β = E [ x 2 ] − α E [ x 1 ] を代入すると
E [ x 2 x 1 ] = α E [ x 1 x 1 ] + E [ x 1 ] ( E [ x 2 ] − α E [ x 1 ] ) = α E [ x 1 2 ] + E [ x 1 ] E [ x 2 ] − α E [ x 1 ] 2 = E [ x 1 ] E [ x 2 ] + α ( E [ x 1 2 ] − E [ x 1 ] 2 ) = E [ x 1 ] E [ x 2 ] + α V [ x 1 ] ⟹ E [ x 2 x 1 ] − E [ x 1 ] E [ x 2 ] = α V [ x 1 ] \begin{align}
E[x_2 x_1]
&= \alpha E[x_1 x_1] + E[x_1] (E[x_2] - \alpha E[x_1])\\
&= \alpha E[x_1^2] + E[x_1] E[x_2] - \alpha E[x_1]^2\\
&= E[x_1] E[x_2] + \alpha (E[x_1^2] - E[x_1]^2)\\
&= E[x_1] E[x_2] + \alpha V[x_1]\\
\implies
E[x_2 x_1] - E[x_1] E[x_2] &= \alpha V[x_1]\\
\end{align} E [ x 2 x 1 ] ⟹ E [ x 2 x 1 ] − E [ x 1 ] E [ x 2 ] = α E [ x 1 x 1 ] + E [ x 1 ] ( E [ x 2 ] − α E [ x 1 ]) = α E [ x 1 2 ] + E [ x 1 ] E [ x 2 ] − α E [ x 1 ] 2 = E [ x 1 ] E [ x 2 ] + α ( E [ x 1 2 ] − E [ x 1 ] 2 ) = E [ x 1 ] E [ x 2 ] + α V [ x 1 ] = α V [ x 1 ] 分散は
V [ x ] = E [ x 2 ] − E [ x ] 2 → E [ x 2 ] = V [ x ] − E [ x ] 2 V[x] = E[x^2] - E[x]^2\\
\to E[x^2] = V[x] - E[x]^2 V [ x ] = E [ x 2 ] − E [ x ] 2 → E [ x 2 ] = V [ x ] − E [ x ] 2 であり、
Σ = E [ ( x − μ ) ( x − μ ) ⊤ ] = E [ x x ⊤ ] − E [ x ] μ ⊤ − μ E [ x ⊤ ] + μ μ ⊤ = E [ x x ⊤ ] − 2 μ μ ⊤ + μ μ ⊤ = E [ x x ⊤ ] − μ μ ⊤ \newcommand{\b}[1]{ \boldsymbol{#1} }
\begin{align}
\b{\Sigma}
&= E[(\b{x} - \b{\mu})(\b{x} - \b{\mu})^\top]\\
&= E[\b{x x}^\top] - E[\b{x}]\b{\mu}^\top - \b{\mu} E[\b{x}^\top] + \b{\mu\mu}^\top\\
&= E[\b{x x}^\top] - 2\b{\mu\mu}^\top + \b{\mu\mu}^\top\\
&= E[\b{x x}^\top] - \b{\mu\mu}^\top\\
\end{align} Σ = E [( x − μ ) ( x − μ ) ⊤ ] = E [ xx ⊤ ] − E [ x ] μ ⊤ − μ E [ x ⊤ ] + μμ ⊤ = E [ xx ⊤ ] − 2 μμ ⊤ + μμ ⊤ = E [ xx ⊤ ] − μμ ⊤ であり、2変数だと
Undefined control sequence: \b at position 15: \begin{align}
\̲b̲{\Sigma}
&=
E\l…
\begin{align}
\b{\Sigma}
&=
E\left[ \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}
\begin{pmatrix} x_1 & x_2 \end{pmatrix} \right]
-
\begin{pmatrix} E[x_1] \\ E[x_2] \end{pmatrix}
\begin{pmatrix} E[x_1] & E[x_2] \end{pmatrix}\\
&=
\begin{pmatrix}
E[x_1^2] & E[x_1 x_2]\\
E[x_2 x_1] & E[x_2^2]
\end{pmatrix}
-
\begin{pmatrix}
E[x_1]^2 & E[x_1]E[x_2]\\
E[x_2]E[x_1] & E[x_2]^2
\end{pmatrix}
\end{align}であることを利用すると
E [ x 2 x 1 ] − E [ x 1 ] E [ x 2 ] = α V [ x 1 ] V [ x 1 x 2 ] = α V [ x 1 ] α = V [ x 1 x 2 ] V [ x 1 ] \begin{align}
E[x_2 x_1] - E[x_1] E[x_2] &= \alpha V[x_1]\\
V[x_1 x_2] &= \alpha V[x_1]\\
\alpha &= \frac{V[x_1 x_2]}{V[x_1]}
\end{align} E [ x 2 x 1 ] − E [ x 1 ] E [ x 2 ] V [ x 1 x 2 ] α = α V [ x 1 ] = α V [ x 1 ] = V [ x 1 ] V [ x 1 x 2 ] 推定では標本の統計量を用いて
α ^ = s 21 s 1 2 β ^ = x ˉ 2 − α ^ x ˉ 1 \begin{align}
\hat{\alpha} &= \frac{s_{21}}{s^2_1}\\
\hat{\beta} &= \bar{x}_2 - \hat{\alpha} \bar{x}_1\\
\end{align} α ^ β ^ = s 1 2 s 21 = x ˉ 2 − α ^ x ˉ 1 とする。
import semopy
import numpy as np
import pandas as pd
# 適当なデータを生成
n = 1000
np.random.seed(0)
x = np.random.uniform(size=n)
e = np.random.normal(size=n)
y = 10 + 3 * x + e
data = pd.DataFrame(dict(y=y, x=x))
# モデルを構築
desc = "y ~ x"
model = semopy.Model(desc)
model.fit(data)
model.inspect()# パス図
semopy.semplot(model, filename="/tmp/path_diagram.png")