Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

回帰分析

単回帰モデル

ある変数x2x_2を別の変数x1x_1から予測する問題を考える。

x2=αx1+β+ex_2 = \alpha x_1 + \beta + e

x1x_1予測変数(predictor variable)、x2x_2基準変数(criterion variable)と呼ばれる確率変数である。

eeは予測誤差を表現する確率変数であり誤差変数といい、以下の性質を仮定する。

  • E[e]=0E[e] = 0(誤差は正負バランスよく出現する)

  • E[ex1]=0E[e x_1] = 0(予測変数の値が大きくても小さくても誤差は大きくなったり小さくなったりする傾向はない)

パラメータの推定には最尤法や最小二乗法やモーメント法などが使われる

モーメント法による単回帰モデルの推定

単回帰モデルの両辺の期待値をとると

E[x2]=αE[x1]+E[β]+E[e]=0    β=E[x2]αE[x1]E[x_2] = \alpha E[x_1] + E[\beta] + \underbrace{ E[e] }_{=0}\\ \implies \beta = E[x_2] - \alpha E[x_1]

となる。

次に、単回帰モデルの両辺にx1x_1をかけてから期待値をとると

E[x2x1]=αE[x1x1]+βE[x1]+E[ex1]=0E[x_2 x_1] = \alpha E[x_1 x_1] + \beta E[x_1] + \underbrace{ E[e x_1] }_{=0}\\

β=E[x2]αE[x1]\beta = E[x_2] - \alpha E[x_1]を代入すると

E[x2x1]=αE[x1x1]+E[x1](E[x2]αE[x1])=αE[x12]+E[x1]E[x2]αE[x1]2=E[x1]E[x2]+α(E[x12]E[x1]2)=E[x1]E[x2]+αV[x1]    E[x2x1]E[x1]E[x2]=αV[x1]\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}

分散は

V[x]=E[x2]E[x]2E[x2]=V[x]E[x]2V[x] = E[x^2] - E[x]^2\\ \to E[x^2] = V[x] - E[x]^2

であり、

Σ=E[(xμ)(xμ)]=E[xx]E[x]μμE[x]+μμ=E[xx]2μμ+μμ=E[xx]μμ\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}

であり、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[x2x1]E[x1]E[x2]=αV[x1]V[x1x2]=αV[x1]α=V[x1x2]V[x1]\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}

推定では標本の統計量を用いて

α^=s21s12β^=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}

とする。

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()
Loading...
# パス図
semopy.semplot(model, filename="/tmp/path_diagram.png")
Loading...