カテゴリー: 科学・工学

ラグランジュの未定乗数法～等式制約条件

/ tau / コメントする

準備

長方形の面積最大化

たとえば2変数の問題として、長方形の周囲長Lを一定として、その面積が最大となる長方形の形状と面積はどのようになるかを考える。この場合、長方形の辺の長さをx, yとすると、問題は以下のように表せる。

(1)   \begin{equation*} \mathrm{maximize} \quad S(x, y) = xy \quad \mathrm{subject~to} \quad x+y = \frac{L}{2} \end{equation*}

これは以下のように代数的に簡単に解けて、答えは正方形とわかる。

(2)   \begin{gather*} S = x \left( \frac{L}{2} - x \right) = -x^2 + \frac{L}{2} x = - \left( x - \frac{L}{4} \right)^2 + \frac{L^2}{16} \\ \max S = \frac{L^2}{16} \quad \mathrm{for} \: x = y = \frac{L}{4} \end{gather*}

ただし変数の数が増えたり、目的関数や制約条件が複雑になると、解析的に解くのが面倒になる。

Lgrangeの未定乗数法による解

解法から先に示す。Lagrangeの未定乗数法では、目的関数L(x, y)に対して以下の問題となる。

(3)   \begin{gather*} \mathrm{maximize} \quad L(x, y, \lambda) = S(x, y) - \lambda g(x, y) = xy - \lambda \left( x + y - \frac{L}{2} \right) \\ \mathrm{subject~to} \quad S(x, y) = xy, \quad g(x, y) = x + y - \frac{L}{2} \end{gather*}

L(x, y, \lambda)を最大化するために、x, y, \lambdaで偏微分した以下の方程式を設定する。

(4)   \begin{gather*} \frac{\partial L}{\partial x} = 0, \quad \frac{\partial L}{\partial y} = 0, \quad \frac{\partial L}{\partial \lambda} = 0 \end{gather*}

これを計算すると

(5)   \begin{gather*} y - \lambda = 0, \quad x - \lambda = 0, \quad x + y - \frac{L}{2}= 0 \\ \therefore \lambda = x = y = \frac{L}{4} \end{gather*}

Lagrangeの未定乗数法の一般形

一般には、変数\boldsymbol{x} = (x_1, \ldots, x_n)について、目的関数f(\boldsymbol{x})を制約条件g(\boldsymbol{x})=0の下で最大化／最小化する問題として与えられる。

(6)   \begin{align*} & \mathrm{maximize/minimize} \quad f(\boldsymbol{x}) \\ & \mathrm{subject~to} \quad g(\boldsymbol{x}) = 0 \end{align*}

この等式制約条件付き最大化／最小化問題は、以下のようにL(\boldsymbol{x}, \lambda)を導入して、連立方程式として表現される。

(7)   \begin{align*} & L(\boldsymbol{x}, \lambda) = f(\boldsymbol{x}) - \lambda g(\boldsymbol{x}) \\ & \frac{\partial L(\boldsymbol{x}, \lambda)}{\partial x_1} = \cdots = \frac{\partial L(\boldsymbol{x}, \lambda)}{\partial x_n} = \frac{\partial L(\boldsymbol{x}, \lambda)}{\partial \lambda} = 0 \end{align*}

例題

例題1：平面と円

平面x + y - 1について、制約条件f(x, y) = x^2 + y^2の下での極値を求める。

(8)   \begin{align*} & \mathrm{minimize} \quad f(x, y) = x + y - 1 \\ & \mathrm{subject~to} \quad x^2 + y^2 = 1 \end{align*}

lagrangeの未定乗数を導入して問題を定式化すると以下のようになる。

(9)   \begin{equation*} L(x, y, \lambda) = x + y - 1 - \lambda (x^2 + y^2 - 1) \end{equation*}

(10)   \begin{align*} &\dfrac{\partial L}{\partial x} = 1 - 2 \lambda x = 0 \\ &\dfrac{\partial L}{\partial y} = 1 - 2 \lambda y = 0 \\ &\dfrac{\partial L}{\partial \lambda} = - x^2 - y^2 + 1 = 0 \end{align*}

この連立方程式を解くと以下のようになり、解として2つの極値を得るが、それらは最大値と最小値に相当する。

(11)   \begin{gather*} x = y = \frac{1}{2 \lambda} \quad \Rightarrow \quad \frac{1}{4 \lambda ^2} + \frac{1}{4 \lambda ^2} = 1 \quad \Rightarrow \quad \lambda = \pm \frac{1}{\sqrt{2}}\\ \therefore \; x = y = \pm \frac{1}{\sqrt{2}} \approx \pm 0.7071\\ \max f(x, y) = \sqrt{2} - 1 \approx 0.414\\ \min f(x, y) = -\sqrt{2} - 1 \approx -2.414 \end{gather*}

これを目的関数のコンターと制約条件の線で表すと以下の通り。

 

例題2：凸関数と直線

下に凸な関数f(x, y) = x^2 + y^2について、直線x + y = 1の制約条件下での最小値を求める。

(12)   \begin{align*} & \mathrm{minimize} \quad f(x, y) = x^2 + y^2 \\ & \mathrm {subject~to} \quad x + y - 1 = 0 \end{align*}

ここでlagrangeの未定乗数を導入して問題を定式化すると以下のようになる。

(13)   \begin{equation*} L(x, y, \lambda) = x^2 + y^2 - \lambda (x + y - 1) \end{equation*}

(14)   \begin{align*} &\dfrac{\partial L}{\partial x} = 2x - \lambda = 0 \\ &\dfrac{\partial L}{\partial y} = 2y - \lambda = 0 \\ &\dfrac{\partial L}{\partial \lambda} = - x - y + 1 = 0 \end{align*}

この連立方程式を解くと以下のようになり、解は最小値1つとなる。

(15)   \begin{gather*} x = y = \frac{\lambda}{2} \quad \Rightarrow \quad \frac{\lambda}{2} + \frac{\lambda}{2} = 1 \quad \Rightarrow \quad \lambda = 1\\ \therefore \; x = y = \frac{1}{2} \quad , \quad \min f(x, y) = \frac{1}{2} \end{gather*}

これを目的関数のコンターと制約条件の線で表すと以下の通り。

なお、これを3次元で表示すると以下のようになる。青い曲面が目的関数で、赤い直線が制約条件となる。最適化問題は、制約条件を満たす曲面上の点(図中、赤い放物線)の最小値を求めることになる。

幾何学的説明

式(7)は以下のように書ける。

(16)   \begin{gather*} \left[ \begin{array}{c} \dfrac{\partial f}{\partial x_1} \\ \vdots \\ \dfrac{\partial f}{\partial x_n} \end{array} \right] = \lambda \left[ \begin{array}{c} \dfrac{\partial g}{\partial x_1} \\ \vdots \\ \dfrac{\partial g}{\partial x_n} \end{array} \right] \\ g(x_1, \ldots, x_n) = 0 \end{gather*}

さらにgradientで表すと

(17)   \begin{gather*} \nabla f = \lambda \nabla g \\ g(x_1, \ldots, x_n) = 0 \end{gather*}

すなわちこの式の解(x_1, \ldots, x_n)は、制約条件であるg(x_1, \ldots, x_n)=0を満足し、その曲線上にある。さらに解の点においてf(x_1, \ldots, x_n)の勾配ベクトルとゼロ平面上におけるg(x_1, \ldots, x_n)の勾配ベクトルが平行になる。これはゼロ平面上の解の点において制約条件の曲線とf(x_1, \ldots, x_n)のコンターのうち特定の曲線が接するのと同義であり、この点は停留点(stationary point)である。

つまりこのような停留点を発見する手順は、制約条件を満たす(制約条件の線上にある)点のうち、その点において目的関数のgradientと制約条件の関数のgradientが平行となる点を求めるということになる。再度これを式で表すと、

(18)   \begin{gather*} \nabla f(\boldsymbol{x}) = \lambda \nabla g(\boldsymbol{x}) \\ g(\boldsymbol{x}) = 0 \end{gather*}

となるが、これをLangrange関数L=f - \lambda gと定義したうえで各変数で偏微分したものをゼロと置いた方程式を解くと表現している。未定乗数λは停留点における目的関数のgradientと制約条件の関数のgradientの比を表している。

λの符号

λの符号に意味があるかどうか。

たとえば、以下の制約条件付き最適化問題を考える。

(19)   \begin{align*} & \mathrm{minimize} \quad f(x, y) = x^2 + y^2 \\ & \mathrm {subject~to} \quad x + y - 2 = 0 \end{align*}

(20)   \begin{gather*} L(x, y, \lambda) = x^2 + y^2 - \lambda(x + y - 2) \\ \frac{\partial L}{\partial x} = \frac{\partial L}{\partial y} = \frac{\partial L}{\partial \lambda} = 0 \\ 2x - \lambda = 2y - \lambda = -x - y + 2 = 0\\ \lambda = 2, \; x=y=1 \end{gather*}

この問題の制約条件を変更して以下のようにした場合。

(21)   \begin{align*} & \mathrm{minimize} \quad f(x, y) = x^2 + y^2 \\ & \mathrm {subject~to} \quad - x - y + 2 = 0 \end{align*}

(22)   \begin{gather*} L(x, y, \lambda) = x^2 + y^2 - \lambda(- x - y + 2) \\ \frac{\partial L}{\partial x} = \frac{\partial L}{\partial y} = \frac{\partial L}{\partial \lambda} = 0 \\ 2x + \lambda = 2y + \lambda = x + y - 2 = 0\\ \lambda = -2, \; x=y=1 \end{gather*}

このように、制約条件の正負を反転するとλの符号が逆になるが解は変わらない。これを表したのが以下の図。

等式制約条件の場合、制約条件の線上でgradientが平行になりさえすればよいので、λ符号(制約条件式の正負)には拘らなくてよい。

停留点が極致とならない例

以下の最適化問題の解をみてみる。

(23)   \begin{align*} & \mathrm{maxmize} \quad f(x, y) = x^3 + y^3 \\ & \mathrm{subject~to} \quad g(x, y) = x - y \end{align*}

(24)   \begin{align*} &L(x, y, \lambda) = x^3 + y^3 - \lambda(x - y) \\ &\left\{ \begin{array}{c} \dfrac{\partial L}{\partial x} = 3x^2 - \lambda = 0\\ \\ \dfrac{\partial L}{\partial y} = 3y^2 + \lambda= 0\\ \\ \dfrac{\partial L}{\partial \lambda} = x - y = 0 \end{array} \right. \end{align*}

第1式から第2式を引き、第3式を適用して、

(25)   \begin{gather*} \3x^2 - 3y^2 - 2\lambda = 0 \\ 3(x + y)(x - y) = 2\lambda = 0 \\ \therefore x = y = \lambda = 0 \end{gather*}

ここでz=f(x, y)の曲面と制約条件g(x, y)=0に対する曲面上の軌跡を描くと下図左のようになる。下図の右はt = x = yとして曲線を表したもので、x = y = 0で勾配は水平になっており、停留点ではあるが極大／極小となっていない。

参考サイト

本記事をまとめるにあたって、下記サイトが大変参考になったことに感謝したい。

ヘッセの標準形～点と平面の距離

/ tau / 2件のコメント

2次元の場合

直線l、点Qの距離を考える。直線と各点の記号、座標を以下のように定義する。

(1)   \begin{equation*} l:ax + by + c = 0 \end{equation*}

Step-1：直線に直交するベクトル

まず、ベクトル(a, b)が直線lに直交することを示す。直線は以下のように媒介変数表示できて、ベクトル(u_x , u_y)は直線に平行なベクトル。

(2)   \begin{equation*} \left\{ \begin{array}{lll} x &=& u_x t + x_0 \\ y &=& u_y t+ y_0 \end{array} \right. \end{equation*}

これを直線の式に代入して、

(3)   \begin{gather*} a (u_x t + x_0) + b (u_y t + y_0) + c = 0 \\ (a u_x + b u_y) t + (a x_0 + b y_0 + c) = 0 \end{gather*}

ここで任意のtに対して上式が成り立つことから、a u_x + b u_y = 0となり、ベクトル(a, b)は直線に垂直であることが示された。

別解：点と直線の距離をパラメーター（媒介変数）によって愚直に求める方法

Step-2：法線ベクトルとの平行条件による導出

与えられた点Pから直線lへの垂線の足をHとすると、\overrightarrow{PH} \parallel (a, b)なので、以下が成り立つ。

(4)   \begin{equation*} \left| \overrightarrow{PH} \cdot (a, b) \right| = \left| \overrightarrow{PH} \right| \cdot | (a, b) | \end{equation*}

ここで(x_h, y_h)は直線l上にあることを考慮し、上式の左辺を以下のように変形できる。

(5)   \begin{equation*} \begin{array}{lll} \left| \overrightarrow{PH} \cdot (a, b) \right| &=& | (x_p - x_h ) a + (y_p - y_h) b | \\ &=& | a x_p + b y_p - (a x_h + b y_h) | \\ &=& | a x_p + b y_p + c | \end{equation*}

これより

(6)   \begin{equation*} \left| \overrightarrow{PH} \right| = \frac{\left| a x_p + b y_p + c \right|}{\sqrt{a^2 + b^2}} \end{equation*}

3次元の場合

3次元平面の式

3次元空間内の平面は、たとえば以下のように表すことができる。

(7)   \begin{equation*} \pi : w_x x + w_y y + w_z z + w_0 = 0 \end{equation*}

一方、3次元平面上の点と法線ベクトル{\boldsymbol n} = (x_x, n_y, n_z)との直行条件から、以下のようにも表現できる。

(8)   \begin{gather*} {\boldsymbol n} ({\boldsymbol x} - {\boldsymbol x_0} ) = 0 \\ (n_x, n_y, n_z) \cdot (x - x_0, y - y_0, z - z_0) = 0 \\ n_x x + n_y y + n_z z + (-n_x x_0 -n_y y_0 -n_z z_0) = 0 \end{gather*}

上記2つの式より、ベクトル{\boldsymbol w} = (w_x, w_y, w_z)は平面に対する法線ベクトルであることがわかる。

法線ベクトルとの平行条件

この法線ベクトルがベクトル\left| \overrightarrow{PH} \right|と平行であることから、

(9)   \begin{equation*} \left| \overrightarrow{PH} \cdot {\boldsymbol w} \right| = \left| \overrightarrow{PH} \right| \cdot \left| {\boldsymbol w} \right| \end{equation*}

上式の左辺は以下のように変形できる。

(10)   \begin{equation*} \begin{array}{lll} \left| \overrightarrow{PH} \cdot {\boldsymbol w} \right| &=& \left| (x_p - x_h, y_p - y_h, z_p - z_h) \cdot (w_x, w_y, w_z) \right| \\ &=& \left| w_x x_p + w_y y_p + w_z z_p - (w_x x_h +w_y y_h + w_z z_h) \right| \\ &=& \left| w_x x_p + w_y y_p + w_z z_p + w_0 \right| \end{array} \end{equation*}

以上のことから、点{\boldsymbol x}_pから三次元平面\piへの距離については、以下で表される。

(11)   \begin{equation*} \left| \overrightarrow{PH} \right| = \frac{\left| w_x x_p + w_y y_p + w_z z_p + w_0 \right|}{| {\boldsymbol w} |} \end{equation*}

多次元の場合

n次元の超平面を以下の式で与える。

(12)   } \begin{equation*} {\boldsymbol w} \cdot {\boldsymbol x} + w_0 = 0 \; \Leftrightarrow \; w_0 + w_1 x_1 + \cdots + w_n x_n = 0 \end{equation*}

このとき、これまでと同様の考え方により、点{\boldsymbol x}_p (x_{p1}, \ldots , x_{pn})と上記の超平面との距離は以下で表される。

(13)   \begin{equation*} \left| \overrightarrow{PH} \right| = \frac{ {\boldsymbol w} \cdot {\boldsymbol x} + w_0 }{ \| {\boldsymbol w} \|} \end{equation*}

 

 

二項分布の平均と分散

/ tau / コメントする

概要

二項分布B(n, p)の平均と分散は以下のようになる。

(1)   \begin{alignat*}{1} E(X) &= np \\ V(X) &= np(1-p) \end{alignat*}

これらを導くのに、有用なテクニックを使っているのでまとめておく。

直接定義式から導く方法

この方法は、平均、分散の定義式から直接導いていく過程で、意図的に二項展開の形に持ち込んでいく方法。

平均

平均の定義から、以下のように変形していく。

(2)   \begin{alignat*}{1} E(X) &= \sum_{k=0}^n k {}_n\mathrm{C}_k p^k (1-p)^{n-k} \\ &= \sum_{k=1}^n k \frac{n!}{(n-k)! k!} p^k (1-p)^{n-k} \\ &= \sum_{k=1}^n n \frac{(n-1)!}{\left( (n-1) - (k-1) \right)! (k-1)!} pp^{k-1} (1-p)^{(n-1)-(k-1)} \end{alignat*}

上式では、各項にkが乗じられていることから、以下の流れで変形している。

  • k=0のとき1番目の項はゼロとなるので、和の開始値をk=1とする
  • 分子にあるkを使って、{}_n C_kの分母においてk! \rightarrow (k-1)!とする
  • これに整合させるため、分母において(n-k)!を((n-1) - (k-1))!と変形
  • さらに組み合わせの式に整合させるため分子をn (n-1)!と変形し、最終的に{}_{n-1} C_{k-1}を引き出している。
  • 後に二項定理を使うため、p, \; 1-pの指数も調整している

ここでk-1 = k'とおくと、カウンターの範囲はk=1 \sim nからk'=0 \sim n-1となることから、

(3)   \begin{alignat*}{1} E(X) &= np \sum_{k'=0}^{n-1} \frac{(n-1)!}{\left( (n-1) - k' \right)! k'!} p^{k'} (1-p)^{(n-1)-k'} \\ &= np (p + 1 - p)^{n-1} \\ &= np \end{alignat*}

上式では、変形した和の部分が二項展開の形になっていることを利用している。

分散

分散については、k^2が各項に乗じられるが、これをk(k-1) + kと変形して、階乗のランクを下げるところがミソ。

(4)   \begin{alignat*}{1} V(X) &= E(X^2) - [E(X)]^2 \\ &= \sum_{k=0}^n k^2 {}_n C_k p^k (1-p)^{n-k} - (np)^2 \\ &= \sum_{k=0}^n \left( k(k-1) + k \right) {}_n C_k p^k (1-p)^{n-k} - (np)^2 \\ &= \sum_{k=2}^n k(k-1) {}_n C_k p^k (1-p)^{n-k} + \sum_{k=1}^n k {}_n C_k p^k (1-p)^{n-k} - (np)^2 \\ \end{alignat*}

上式で、1項目はk(k-1)が乗じられているのでカウンターをk=2から、2項目は同じくk=1からとしている。

ここで1項目についてk'' = k - 2と置いて、平均の時と同じ考え方で以下のように変形。

(5)   \begin{alignat*}{1} &\sum_{k=2}^n k(k-1) {}_n C_k p^k (1-p)^{n-k} \\ &= \sum_{k=2}^n k(k-1) \frac{n!}{(n-k)! k!} p^k (1-p)^{n-k} \\ &= \sum_{k=2}^n n(n-1) \frac{(n-2)!}{((n-2) -(k-2))! (k-2)!} p^2 p^{k-2} (1-p)^{(n-2)-(k-2)} \\ &= n(n-1) p^2 \sum_{k'=0}^{n-2} \frac{(n-2)!}{((n-2) -k'')! k''!} p^{k''} (1-p)^{(n-2)-k''} \\ &= n(n-1) p^2 ( p + (1-p))^{n-2} \\ &= n(n-1) p^2 \end{alignat*}

2項目については、k' = k-1とおいて、

(6)   \begin{alignat*}{1} &\sum_{k=1}^n k {}_n C_k p^k (1-p)^{n-k} \\ &= \sum_{k=1}^n n \frac{(n-1)!}{\left( (n-1) - (k-1) \right) ! (k-1)!} p p^{k-1} (1-p)^{(n-1)-(k-1)} \\ &= np \sum_{k'=0}^{n-1} \frac{(n-1)!}{\left( (n-1) - k'\right) ! k'!} p p^{k'} (1-p)^{(n-1)-k'} \\ &= np \left( (p + (1-p) \right) ^{n-1} \\ &= np \end{alignat*}

以上を併せて、

(7)   \begin{alignat*}{1} V(X) &= n(n-1) p^2 + np - (np)^2 \\ &= np(1-p) \end{alignat*}

微分による方法

この方法は、kp^k, \; k^2 p^kの形に着目して、全事象の式を微分する方法。式展開が素直であり、平均の式を微分した結果がそのまま分散の式になってしまうところが美しい。

平均

二項分布の全確率の和は1となる。

(8)   \begin{equation*} \sum_{k=0}^n {}_n \mathrm{C}_k p^k (1-p)^{n-k} = (p + 1 - p)^n = 1 \end{equation*}

この式の両辺をpで微分する。

(9)   \begin{gather*} \sum_{k=0}^n {}_n \mathrm{C}_k \left( k p^{k-1} (1-p)^{n-k} - (n-k)p^k (1-p)^{n-k-1} \right) = 0 \\ \sum_{k=0}^n {}_n \mathrm{C}_k p^{k-1} (1-p)^{n-k-1} \left( k(1-p) - (n-k)p \right) = 0 \\ \sum_{k=0}^n {}_n \mathrm{C}_k p^{k-1} (1-p)^{n-k-1} (k - np) = 0 \\ \end{gather*}

両辺にp(1-p)をかける。

(10)   \begin{gather*} \sum_{k=0}^n {}_n \mathrm{C}_k p^k (1-p)^{n-k} (k - np) = 0 \\ \sum_{k=0}^n k {}_n \mathrm{C}_k p^k (1-p)^{n-k} = np \sum_{k=0}^n {}_n \mathrm{C}_k p^k (1-p)^{n-k}\\ \therefore E(X) = np \end{gather*}

分散

式(10)をもう一度pで微分する。

(11)   \begin{gather*} \sum_{k=0}^n k {}_n \mathrm{C}_k p^k (1-p)^{n-k} = np \\ \sum_{k=0}^n k {}_n \mathrm{C}_k \left( kp^{k-1} (1-p)^{n-k} -(n-k)p^k (1-p)^{n-k-1} \right) = n\\ \sum_{k=0}^n k {}_n \mathrm{C}_k p^{k-1} (1-p)^{n-k-1} \left( k (1-p) -(n-k)p \right) = n\\ \sum_{k=0}^n k {}_n \mathrm{C}_k p^{k-1} (1-p)^{n-k-1} (k - np) = n\\ \end{gather*}

両辺にp(1-p)をかける。

(12)   \begin{gather*} \sum_{k=0}^n k {}_n \mathrm{C}_k p^k (1-p)^{n-k} (k - np) = np(1-p) \\ \sum_{k=0}^n k^2 {}_n \mathrm{C}_k p^k (1-p)^{n-k} - np \sum_{k=0}^n k {}_n \mathrm{C}_k p^k (1-p)^{n-k}= np(1-p) \\ E(X^2) - (np)^2 = np(1-p) \\ \therefore V(X) = np(1-p) \end{gather*}

 

回帰分析～重回帰

/ tau / コメントする

概要

単回帰が(X, Y)への線形関係を扱うのに対して、重回帰は複数の説明変数X_k \; (k=1 \ldots m)とYの線形関係を扱う。

データ(X_1, \ldots , X_m, Y) = (x_{1i}, \ldots, x_{mi}, y_i), (i = 1 \ldots n)に対して以下の線形式で最も説明性の高いものを求める。

(1)   \begin{equation*} \hat{y} = w_0 + w_1 x_1 + \cdots + w_m x_m \end{equation*}

そのために、データとその推測値の残差\hat{y_i} - y_iの平方和が最小となるような係数w_0, w_1, \ldots, w_mを最小二乗法により求める。

定式化

説明変数とターゲットのデータセットn組が次のように得られているとする。

(2)   \begin{equation*} x_{1i}, \ldots , x_{mi}, y_i \quad {\rm where} \quad i = 1, \ldots , n \end{equation*}

残差の平方和については

(3)   \begin{equation*} \min \sum_{i=1}^n (\hat{y} -y_i)^2 \quad {\rm where} \quad \hat{y}_i = w_0 + w_1 x_{1i} + \cdots w_m x_{mi} \end{equation*}

ここで残差を最小化するw_0, w_1, \ldots, w_mを求めるために、それぞれで偏微分する。

(4)   \begin{equation*} \begin{array}{c} \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial w_0} \sum_{i=1}^n (w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i)^2 = 0 \\ \displaystyle \frac{\partial}{\partial w_1} \sum_{i=1}^n (w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i)^2 = 0 \\ \vdots \\ \displaystyle \frac{\partial}{\partial w_m} \sum_{i=1}^n (w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i)^2 = 0 \\ \end{array} \right. \end{equation*}

変形すると、

(5)   \begin{equation*} \left\{ \begin{array}{l} \displaystyle \sum_{i=1}^n 2(w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i) = 0 \\ \displaystyle \sum_{i=1}^n 2(w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i) x_{1i} = 0 \\ \vdots \\ \displaystyle \sum_{i=1}^n 2(w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i) x_{mi} = 0 \\\end{array} \right. \end{array} \end{equation*}

更にこれを展開するのに、S_k = \sum_{i=1}^n x_{ki}, S_{kl} = \sum_{i=1}^n x_{ki} x_{li}, S_y = \sum_{i=1}^n y_i, S_{ky} = \sum_{i=1}^n x_{ki} y_iとおいて

(6)   \begin{equation*} \left\{ \begin{array}{l} n w_0 + w_1 S_1+ \cdots + w_m S_m - S_y = 0 \\ w_0S_1 + w_1 S_{11} + \cdots + w_m S_{1m} - S_{1y} = 0 \\ \vdots \\ w_0 S_m + w_1 S_{m1} + \cdots + w_m S_{mm} - S_{my} = 0 \\ \end{array} \right. \end{array} \end{equation*}

これを行列形式で表示すると

(7)   \begin{equation*} \left( \begin{array}{cccc} n & S_1 & \cdots & S_m \\ S_1 & S_{11} & \cdots & S_{1m} \\ \vdots & \vdots & & \vdots \\ S_m & S_{m1} & \cdots & S_{mm} \end{array} \right) \left( \begin{array}{c} w_0 \\ w_1 \\ \vdots \\ w_m \end{array} \right) = \left( \begin{array}{c} S_y \\ S_{1y} \\ \vdots \\ S_{my} \end{array} \right) \end{equation*}

解の導出

式(7)の連立方程式を解けば各係数を得ることができる。ここで左辺の係数行列は以下のようにも表せる。

(8)   \begin{equation*} \left( \begin{array}{cccc} 1 & 1 & \cdots & 1 \\ x_{11} & x_{12} & \cdots & x_{1n} \\ \vdots & \vdots & & \vdots \\ x_{m1} & x_{m2} & \cdots & x_{mn} \end{array} \right) \left( \begin{array}{cccc} 1 & x_{11} & \cdots & x_{m1} \\ 1 & x_{12} & \cdots & x_{m2} \\ \vdots & \vdots & & \vdots \\ 1 & x_{1n} & \cdots & x_{mn} \end{array} \right) = \boldsymbol{X}^T \boldsymbol{X} \end{equation*}

また右辺は次のように表せる。

(9)   \begin{equation*} \left( \begin{array}{cccc} 1 & 1 & \cdots & 1 \\ x_{11} & x_{12} & \cdots & x_{1n} \\ \vdots & \vdots & & \vdots \\ x_{m1} & x_{m2} & \cdots & x_{mn} \end{array} \right) \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array} \right) = \boldsymbol{X}^T \boldsymbol{y} \end{equation*}

したがって式(7)は以下のように表せる。

(10)   \begin{equation*} \boldsymbol{X}^T \boldsymbol{X} \boldsymbol{w} = \boldsymbol{X}^T \boldsymbol{y} \end{equation*}

これを解くと

(11)   \begin{equation*} \boldsymbol{w} = ( \boldsymbol{X}^T \boldsymbol{X} )^{-1} \boldsymbol{X}^T \boldsymbol{y} \end{equation*}

ベクトル・行列による表現と解

定式化の段階からベクトル・行列を用いてみる。式(1)の各変数を全てベクトル形式で表すと以下のようになる。

(12)   \begin{equation*} \hat{\boldsymbol{y}} = \boldsymbol{x}^T \boldsymbol{w} \end{equation*}

ここで、

(13)   \begin{equation*} \boldsymbol{x} = \left( \begin{array}{c} 1 & x_1 & \vdots & x_n \end{array} \right) \quad , \quad \boldsymbol{w} = \left( \begin{array}{c} w_0 & w_1 & \vdots & w_n \end{array} \right) \end{equation*}

複数データからなるデータセットに対して、式(8)の行列表現を用いると、以下のように表現できる。

(14)   \begin{equation*} \hat{\boldsymbol{y}} = \boldsymbol{X} \boldsymbol{w} \end{equation*}

各データの残差については以下のようになる。

(15)   \begin{equation*} \hat{\boldsymbol{y}} - \boldsymbol{y} = \boldsymbol{X} \boldsymbol{w} - \boldsymbol{y} \end{equation*}

ここで残差の二乗和を損失関数Lとして、これを行列で表現すると以下のようになる。

(16)   \begin{align*} L &= \left( \hat{\boldsymbol{y}} - \boldsymbol{y} \right)^T \left( \hat{\boldsymbol{y}} - \boldsymbol{y} \right) \\ &= \left( \boldsymbol{Xw} - \boldsymbol{y} \right)^T \left( \boldsymbol{Xw} - \boldsymbol{y} \right) \\ &= \left( \boldsymbol{w}^T \boldsymbol{X}^T - \boldsymbol{y}^T \right) \left( \boldsymbol{Xw} - \boldsymbol{y} \right) \\ &= \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{Xw} - \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{y} - \boldsymbol{y}^T \boldsymbol{Xw} + \boldsymbol{y}^T \boldsymbol{y} \end{align*}

上式の第2項、第3項については以下のように整理できる。

(17)   \begin{align*} - \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{y} - \boldsymbol{y}^T \boldsymbol{Xw} = - \left( \boldsymbol{Xw} \right)^T \boldsymbol{y} - \boldsymbol{y}^T \boldsymbol{Xw} = -2 \boldsymbol{y}^T \boldsymbol{Xw} \end{align*}

損失関数Lはスカラーであり、これを最小にするためにベクトルwで微分し、それらの値がゼロとなるような方程式とする。

(18)   \begin{align*} \frac{dL}{d\boldsymbol{w}} = \frac{d}{d\boldsymbol{w}} \left( \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{Xw} - 2 \boldsymbol{y}^T \boldsymbol{Xw} + \boldsymbol{y}^T \boldsymbol{y} \right) = 2 \boldsymbol{X}^T \boldsymbol{Xw} -2 \boldsymbol{X}^T \boldsymbol{y} = \boldsymbol{0} \end{align*}

なお第1項については、式(7)と式(8)より\boldsymbol{X}^T\boldsymbol{X}は対象行列なので転置しても同じ行列となり、行列の微分の公式から以下のようになることを利用している。

(19)   \begin{align*} \frac{d}{d\boldsymbol{x}} \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{Xw} = \left( \boldsymbol{X}^T \boldsymbol{X} + \left(\boldsymbol{X}^T \boldsymbol{X} \right)^T \right) \boldsymbol{w} = 2 \boldsymbol{X}^T \boldsymbol{X} \boldsymbol{w} \end{align*}

式(18)をwについて解いて以下を得る。

(20)   \begin{equation*} \boldsymbol{w} = \left( \boldsymbol{X}^T \boldsymbol{X} \right)^{-1} \boldsymbol{X}^T \boldsymbol{y} \end{equation*}

多重共線性

式(7)からw_0を消去する。

(21)   \begin{equation*} \left( \begin{array}{ccc} S_{11} - \dfrac{{S_1}^2}{n} & \cdots & S_{1m} - \dfrac{S_1 S_m}{n}\\ \vdots & & \vdots \\ S_{m1} - \dfrac{S_m S_1}{n} & \cdots & S_{mm} - \dfrac{{S_m_}^2}{n} \end{array} \right) \left( \begin{array}{c} w_1 \\ \vdots \\ w_m \end{array} \right) = \left( \begin{array}{c} S_{1y} - \dfrac{S_1 S_y}{n} \\ \vdots \\ S_{my} - \dfrac{S_m S_y}{n} \end{array} \right) \end{equation*}

これを分散・共分散で表すと以下の通り。ただし以下でV_{ii}=\rm{Var}(X_i)、V_{ij} = \rm{Cov}(X_i, X_j)と表している。

(22)   \begin{equation*} \left( \begin{array}{ccc} V_{11} & \cdots & V_{1m}\\ \vdots & & \vdots \\ V_{m1} & \cdots & V_{mm} \end{array} \right) \left( \begin{array}{c} w_1 \\ \vdots \\ w_m \end{array} \right) = \left( \begin{array}{c} V_{1y} \\ \vdots \\ V_{my} \end{array} \right) \end{equation*}

ここでx_j = a x_i + bと2つの間に完全な線形関係がある場合、分散・共分散の性質から以下の関係が成り立つ。

(23)   \begin{equation*} V_{jj} = a^2V_{ii}, \; V_{ji} = V_{ij} = aV_{ii}, \; V_{jk} = V_{kj} = aV_{ji} = aV_{ij} \end{equation*}

これらを式(22)の係数行列に適用すると、

(24)   \begin{align*} &\left[ \begin{array}{ccccccc} V_{11} & \cdots & V_{1i} & \cdots & V_{1j} & \cdots & V_{1m}\\ \vdots && \vdots && \vdots && \vdots\\ V_{i1} & \cdots & V_{ii} & \cdots & V_{ij} & \cdots & V_{im}\\ \vdots && \vdots && \vdots && \vdots\\ V_{j1} & \cdots & V_{ji} & \cdots & V_{jj} & \cdots & V_{jm}\\ \vdots && \vdots && \vdots && \vdots\\ V_{m1} & \cdots & V_{mi} & \cdots & V_{mj} & \cdots & V_{mm} \end{array} \right]\\ &= \left[ \begin{array}{ccccccc} V_{11} & \cdots & V_{1i} & \cdots & aV_{1i} & \cdots & V_{1m}\\ \vdots && \vdots && \vdots && \vdots\\ V_{i1} & \cdots & V_{ii} & \cdots & aV_{ii} & \cdots & V_{im}\\ \vdots && \vdots && \vdots && \vdots\\ aV_{i1} & \cdots & aV_{ii} & \cdots & a^2V_{ii} & \cdots & aV_{im}\\ \vdots && \vdots && \vdots && \vdots\\ V_{m1} & \cdots & V_{mi} & \cdots & aV_{mi} & \cdots & V_{mm} \end{array} \right] \end{align*}

係数行列の行(列)が一次従属となっていて、その行列式がゼロとなることから、連立方程式は解を持たない。変数間の関係が線形でなくても相関がかなり高い場合は行列式がゼロに近づき、解が不安定になる。

numpy – 行列（ndarray)

/ tau / コメントする

ベクトルと行列の定義

リテラル

ベクトルはnp.array()で引数にリストを指定して定義。

行列は同じくnp.array()で引数に二次元配列のリストを指定して定義。

単位行列

numpy.identity(n)でn×nの単位行列を生成。

転置

行列の転置にはtranspose()メソッドを使う。代替として.Tとしてもよい。

一次元配列で定義したベクトルにはtranspose()は効かない。列ベクトルに変換するにはreshape()メソッドを使う（reshape(行数, 列数)）。

演算

定数倍

ベクトル・行列の定数倍は、各要素の定数倍。

加減

同じ要素数のベクトル、同じ次元の行列同士の下限は要素同士の加減

ベクトルの内積

同じ要素数のベクトルの内積（ドット積）はnp.dot()で計算。

    \begin{equation*} {\bf a} \cdot {\bf b} = \sum_{i=1}^n a_i b_i \end{equation*}

*演算子を使うと、要素ごとの積になる。

行列の積

行列同士の積もnp.dot()で計算。l×m行列とm×n行列の積はl×n行列になる。

    \begin{equation*} \left( \begin{array}{ccc} a_{11} & \cdots & a_{1m} \\ \vdots & a_{ij} & \vdots \\ a_{l1} & \cdots & a_{lm} \\ \end{array} \right) \cdot \left( \begin{array}{ccc} b_{11} & \cdots & b_{1n} \\ \vdots & b_{jk} & \vdots \\ b_{m1} & \cdots & b_{mn} \\ \end{array} \right) = \left( \sum_{j=1}^m a_{ij} b_{jk} \right) \end{equation*}

次元数が整合しないとエラーになる。

行ベクトルと行列の積は、ベクトルを前からかけてok。

    \begin{equation*} (1, 2, 3) \left( \begin{array}{ccc} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{array} \right) = (30, 36, 42) \end{equation*}

行列と列ベクトルの積は、一次元配列のベクトルをreshape()で列ベクトルに変換してから。

    \begin{equation*} \left( \begin{array}{ccc} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{array} \right) \left( \begin{array}{c} 1 \\ 2 \\ 3 \end{array} \right) = \left( \begin{array}{c} 14 \\ 32 \\ 50 \end{array} \right) \end{equation*}

なお、np.dot()の代わりに演算子@が使える。ベクトル同士なら内積、少なくともいずれか一つが行列なら行列積。

numpy.linalgパッケージ

行列式

行列式はnumpy.linalgパッケージのdet()関数で得られる。linalgは”linear algebra”の略で、慣例としてLAという名前で代替される。

逆行列

逆行列はnumpy.linalgパッケージのinv()関数で得られる。

固有値・固有ベクトル

正方行列の固有値と固有ベクトルを、eig()関数で得ることができる(行列の固有値・固有ベクトルの例題を用いた)。

結果は固有値が並んだベクトルと固有ベクトルが並んだ配列で、それぞれを取り出して利用する。なお、固有ベクトルはノルムが1となるように正規化されている。

注意が必要なのは固有ベクトルの方で、各固有ベクトルは配列の列ベクトルとして並んでいる。固有ベクトルを取り出す方法は2通り。

固有値に対応するサフィックスで列ベクトルを取り出す。この方法はnumpyの公式ドキュメントにも以下のように書かれている。

v(…, M, M) array
The normalized (unit “length”) eigenvectors, such that the column v[:,i] is the eigenvector corresponding to the eigenvalue w[i].

固有ベクトルの配列を転置して、行ベクトルの並びにする。

 

行列の固有値・固有ベクトル

/ tau / コメントする

概要

行列{\rm A}の固有値・固有ベクトルは以下で定義される。

(1)   \begin{equation*} \boldsymbol{Ax} = \lambda \boldsymbol{x} \end{equation*}

これを以下のように変形する。

(2)   \begin{equation*} (\boldsymbol{A} - \lambda \boldsymbol{I} ) \boldsymbol{x} = {\bf 0} \end{equation*}

この方程式が解をもつためには、以下の条件が必要。

(3)   \begin{equation*} | \boldsymbol{A} - \lambda \boldsymbol{I} | = 0 \end{equation*}

例題

以下の行列に対する固有値、固有ベクトルを求める。

(4)   \begin{equation*} \boldsymbol{A} = \left( \begin{array}{cc} 3 & 1 \\ 2 & 4 \end{array} \right) \end{equation*}

この行列に対する固有値方程式は以下の通り。

(5)   \begin{equation*} | \boldsymbol{A} - \lambda \boldsymbol{I} | = \left| \left( \begin{array}{cc} 3 & 1 \\ 2 & 4 \end{array} \right) - \lambda \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right) \right| = \left| \begin{array}{cc} 3 - \lambda & 1 \\ 2 & 4 - \lambda \end{array} \right| = 0 \end{equation*}

これを解くと、

(6)   \begin{align*} & (3 - \lambda) (4 - \lambda) - 2 = 0 \\ & \lambda ^2 - 7 \lambda + 10 = 0 \\ & (\lambda - 2)(\lambda - 5) = 0 \\ & \lambda = 2, \; 5 \end{align*}

次に、各固有値に対する固有ベクトルを求める。

まず\lambda = 2に対しては、

(7)   \begin{gather*} \left( \begin{array}{cc} 3 - 2 & 1 \\ 2 & 4 - 2 \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{cc} 1 & 1 \\ 2 & 2 \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{c} 0 \\ 0 \end{array} \right) \\ \Rightarrow \; y = -x \\ \therefore \; \boldsymbol{x} = (t, -t) \end{gather*}

確認してみると、

(8)   \begin{gather*} \boldsymbol{Ax} = \left( \begin{array}{cc} 3 & 1 \\ 2 & 4 \end{array} \right) \left( \begin{array}{c} t \\ -t \end{array} \right) = \left( \begin{array}{c} 2t \\ -2t \end{array} \right) , \quad \lambda \boldsymbol{x} = 2 \left( \begin{array}{c} t \\ -t \end{array} \right) = \left( \begin{array}{c} 2t \\ -2t \end{array} \right) \end{gather*}

また\lambda = 5に対しては、

(9)   \begin{gather*} \left( \begin{array}{cc} 3 - 5 & 1 \\ 2 & 4 - 5 \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{cc} -2 & 1 \\ 2 & -1 \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{c} 0 \\ 0 \end{array} \right) \\ \Rightarrow \; y = 2x \\ \therefore \; \boldsymbol{x} = (t, 2t) \end{gather*}

こちらも確認してみると、

(10)   \begin{gather*} \boldsymbol{Ax} = \left( \begin{array}{cc} 3 & 1 \\ 2 & 4 \end{array} \right) \left( \begin{array}{c} t \\ 2t \end{array} \right) = \left( \begin{array}{c} 5t \\ 10t \end{array} \right) , \quad \lambda \boldsymbol{x} = 5 \left( \begin{array}{c} t \\ 2t \end{array} \right) = \left( \begin{array}{c} 5t \\ 10t \end{array} \right) \end{gather*}

なお、固有ベクトルを数値で表現する際、ノルムが1となるように正規化することが多い。

(11)   \begin{gather*} \boldsymbol{x} \rightarrow \frac{\boldsymbol{x}}{\| {\boldsymbol{x}} \|} \end{gather*}

上の例で固有値ベクトルを正規化すると以下の通り。

(12)   \begin{gather*} \frac{(t, -t)}{\sqrt{t^2 + t^2}} = \frac{(1, -1)}{\sqrt{2}} \approx (0.7071, -0.7071) \\ \frac{(t, 2t)}{\sqrt{t^2 + 4t^2}} = \frac{(1, 2)}{\sqrt{5}} \approx (0.4472, 0.8944) \end{gather*}

 

単回帰分析～コンビニエンスストア

/ tau / コメントする

概要

Pythonを使った統計計算と図示の練習のため、コンビニエンスストアで単回帰分析をやってみた。

コンビニの店舗数は「商業動態統計年報」の2016年データを使い、説明変数として2015年の国勢調査人口、2018年の国土地理院による都道府県面積、2017年道路統計年報の道路実延長データを使った。

計算コードは以下の通り。

人口との関係

コンビニ店舗数と人口の散布図と回帰式を以下に示す。

相関係数が極めて高いのは当然で、やはりコンビニの出店計画には人口ファクターが強く作用していることがわかる。

定数項bが店舗数のオーダーに比べてほぼゼロというのも興味深い。

係数aの値からは5万人弱で1店舗ということになるが、人口10万くらいの都市で2店舗しかないことになり、ちょっと少ないような気がする。もしかすると、一定規模以上の市町村単位や都市単位くらいで層化して出店計画を立てているのかもしれない。

面積との関係

店舗数と面積の関係

次に店舗数と都道府県面積の関係を見てみた。

この結果はかなりはずれで、データを見ても人口が少なく面積が群を抜いて大きい北海道の影響を大きく受けている。

ここで、面積が極端に大きい北海道（83423.83㎢～1位、2906店舗～5位）と面積が小さいが集積度が極端に高い東京都（2193.96㎢～45位、7003店舗～1位）の2つを除いて計算してみる。

これはさらにおかしな結果で、面積が小さいほど店舗数が多いことになる。

考えてみれば、集客数を期待するなら人口が集積している地域が有利だから、人口密度に比例する可能性を考えた方がいいのかもしれない。もし面積が小さい県の方が集積度が高いと想定すると、面積だけを取り出したときに逆の関係になるとも考えられるが、相関係数や決定係数が小さすぎるので考察は難しい。

人口密度

以下は人口密度との関係。

今度はかなりきれいに相関の高さが出ている。

直接的な計算式に入れているかどうかわからないが、GISなどで出店計画を立てるとしたら、人口密度の高いエリアを選んでいくだろうことが想定される。

ただ、店舗数は人口などといった売り上げに直結するデータから導かれるのが普通で、人口密度が高くても人口が少なければ出店インセンティブにはならない。

人口と人口密度の関係

試しに人口を説明変数、人口密度を被説明変数として両者の関係を見てみると、驚くことに「人口が多いほど人口密度が高くなる（あるいはその逆）」という関係になる。

ここから先は人口論や地域論になりそうなので置いておくが、少なくとも日本においては、「狭いところほど人が集まっている傾向がある」ということになりそうである。

もちろんこれは他の国でも一般に当てはまることかもしれないが、朝のラッシュ時に特定の車両に無理やり乗り込んでいる割に離れた車両がすいているとか、1本電車を遅らせたらガラガラだったとか、そのあたりの行動パターンを見ていると、何となく日本に特有のような気がする。

道路延長との関係

最後に道路延長との関係を見てみる。

東京のように稠密な都市は例外とすると、概ね関係はありそうである。ただし相関係数、決定係数は高くない。

コンビニ店舗が道路の利便性に依っていることは推測できるが、やはり人口という売り上げ直結のデータに比べると関係は弱い。

高速道路の延長についても見てみたが、こちらはほとんど関係は見られなかった。

ただ、高速道路の延伸に伴ってコンビニエンスストアの店舗数が伸びているようであり、マクロな延長というよりも物流上のインパクトが大きいことは予想される。

 

相関係数

/ tau / コメントする

相関係数の定義

相関係数は、多数のデータの組がどの程度線形に近い性質を持つかを表す値で、以下で定義される。

(1)   \begin{eqnarray*} r &=& \frac{{\rm Cov}(X, Y)}{\sqrt{V(X) \cdot V(Y)}} \\ &=& \frac{E(X - \overline{X})(Y - \overline{Y})} {\sqrt{ ( E\left[ (X - \overline{X})^2 \right] E\left[ (Y - \overline{Y})^2 \right] }} \\ &=& \frac{\displaystyle \sum_{i=1}^{n}(x_i - \overline{x})(y_i - \overline{y})} {\left[ \displaystyle \sum_{i=1}^{n} (x_i - \overline{x})^2 \displaystyle \sum_{i=1}^{n} (y_i - \overline{y})^2 \right]^{1/2}} \end{eqnarray*}

相関係数の線形変換

変数を線形変換した場合の相関係数は分散・共分散の性質から、以下のようになり、元のままか符号が反転する。

(2)   \begin{eqnarray*} r' &=& \frac{{\rm Cov}(aX+b, cY+d)}{\sqrt{V(aX+b) \cdot V(cY+d)}} \\ &=& \frac{ac{\rm Cov}(X, Y)}{\sqrt{a^2 V(X) \cdot c^2 V(Y)}} \\ &=& \frac{ac}{|ac|} \cdot r \end{eqnarray*}

完全線形関係の場合の相関係数

XとYが完全な線形関係にある場合、相関係数は1または-1になる。このとき、傾きの大きさや、平行移動は相関係数に影響しない。

(3)   \begin{eqnarray*} r &=& \frac{ {\rm Cov}(X, aX+b) } { \left[ V(X) \cdot V(aX+b) \right]^{1/2} } \\ &=& \frac{ a {\rm Cov}(X, X) } { \left[ V(X) \cdot a^2 V(X) \right]^{1/2} } \\ &=& \frac{a}{\sqrt{a^2}} \cdot \frac{ {\rm Cov}(X, X) } { \left[ V(X) \cdot V(X) \right]^{1/2} } \\ &=& \frac{a}{|a|} \cdot \frac{ V(X) } { V(X) } \\ &=& 1 \; {\rm or} \; -1 \end{eqnarray*}

いろいろな分布の相関係数

完全な線形関係

以下のコードで確認。

相関係数は傾きや平行移動に対して影響を受けず、増加関数なら1、減少関数なら－1になることがわかる。

線形性が強い関係

線形関数の値に対して、乱数でばらつきを与えた場合の相関係数の違いを示す。ばらつきが大きい方が相関係数は小さくなる。

負の線形性が強い場合は、相関係数がマイナスになる。

放物線（強い関係があるのに相関が低いケース）

以下のような放物線では、XとYにきちんとした数学的関係があるのに、相関係数がゼロに近くなる。

相関係数はXとYが単調増加／単調減少の度合いが強いほど、また線形関係に近いほど1に近くなるが、それ以外の関係が強い場合にはそれを補足できない場合がある。

反比例（負の線形性に見えてしまう場合）

以下は反比例関数の場合。

関数の形状や範囲に寄るが、この場合は相関係数の絶対値が0.8以上と1に近く、これだけ見ると負の線形性が強そうに見える。

対数関数（正の線形性に見えてしまう場合）

対数関数の場合。この場合は0.9以上とかなり強い線形性を示唆している。

相関係数に関する注意

本来の関係との乖離

先にみたように、線形関係ではないが数学的な関係を持つ場合に、相関係数からは全く関係がない、元の関係とは異なり線形関係を持つ、といった解釈になることがある。

相関係数が高い場合に、線形回帰式などで物事を予測する際には注意が必要。

変数間に解析的な関係が見いだせるならそれを重視すべきであり、よしんばそれがわからないにしても、定義域の範囲で「ある程度は当たる」程度に考えておくべきか。

因果関係

堂々と間違えられるケースが、「科学的な」記事やマスメディアなどでよくみられる。

相関係数は「変数間の単調な増加／減少傾向が強いかどうか」だけを示すもので、必ずしも因果関係を示唆しない。

気温が高いとビールはよく売れるがおでんは売れない。その二つに負の相関があるからといって、「ビール好きはおでんが嫌い」と言えないが、形を変えてこのような解釈がなされる恐れがある。

もともとのメカニズムで因果関係が示唆されていて、その上で相関係数の大きさを論じるなら意義もあるが、その場合でも、事象に対する寄与度などをよく考えておかないと「それだけが原因」と考えるような間違いを犯すことになる。

 

回帰分析～単回帰

/ tau / コメントする

概要

(X,\: Y) = (x_i,\: y_i) \; (i=1 \ldots n)のデータからなる複数の標本に対して、線形式で最も説明性の高いものを求める。

そのために、各データのx_iに線形式を適用して得られた値\hat{y_i} = a x_i + bとy_iの差（残差）の平方和が最小となるような係数a,\: bを求める（最小二乗法）。

定式化

(1)   \begin{equation*} \sum_{i=1}^n (\hat{y_i} - y_i)^2 \; \rightarrow \; \min \quad {\rm where} \; \hat{y_i} = a x_i + b \end{equation*}

ここで残差を最小化するa,\: bを求めるために、それぞれで偏微分する。

(2)   \begin{equation*} \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial a} \sum_{i=1}^n (a x_i + b - y_i)^2 = 0 \\ \displaystyle \frac{\partial}{\partial b} \sum_{i=1}^n (a x_i + b - y_i)^2 = 0 \end{array} \right. \quad \Rightarrow \quad \left\{ \begin{array}{l} \displaystyle \sum_{i=1}^n 2(a x_i + b - y_i)x_i = 0 \\ \displaystyle \sum_{i=1}^n 2(a x_i + b - y_i) = 0 \end{array} \right. \end{equation*}

展開して

(3)   \begin{equation*} \left\{ \begin{array}{l} \displaystyle a \sum_{i=1}^n x_i^2 + b \sum_{i=1}^n x_i - \sum_{i=1}^n x_i y_i = 0 \\ \displaystyle a \sum_{i=1}^n x_i + bn - \sum_{i=1}^n y_i = 0 \end{array} \right. \end{equation*}

各和の記号を新たに定義し、行列表示で表すと、

(4)   \begin{equation*} \left[ \begin{array}{cc} S_{xx} & S_x \\ S_x & n \end{array} \right] \left[ \begin{array}{c} a \\ b \end{array} \right] = \left[ \begin{array}{c} S_{xy} \\ S_y \end{array} \right] \end{equation*}

回帰式の係数

式(4)を解くと以下を得る。

(5)   \begin{equation*} \left[ \begin{array}{c} a \\ b \end{array} \right] = \left[ \begin{array}{c} \dfrac{n S_{xy} - S_x S_y}{n S_{xx} - {S_x}^2} \\ \dfrac{S_{xx} S_y - S_x S_{xy}}{n S_{xx} - {S_x}^2} \end{array} \right] \end{equation*}

あるいは式(5)については以下のようにも表せる。ただし以下の式で、\sigma_{XY} = {\rm Cov}(X, Y)、{\sigma_X}^2 = {\rm Var}(X)

(6)   \begin{equation*} a = \frac{\displaystyle \sum_{i=1}^n (x_i - \overline{X})(y_i - \overline{Y})}{\displaystyle \sum_{i=1}^n (x_i - \overline{X})^2} =\frac{\sigma_{XY}}{{\sigma _X}^2} \end{equation*}

(7)   \begin{equation*} b = \overline{Y} - a \overline{X} \end{equation*}

決定係数

全変動・回帰変動・残差変動

x_iの回帰式による推定値を\hat{y}_i = a x_i + bで表す。

ここで、y_i,\: \hat{y},\: \overline{Y}によって、以下の変動が定義される。

y_i - \overline{Y}
全変動。平均からのデータのばらつき。
\hat{y}_i - \overline{Y}
回帰変動。平均からのばらつきのうち、回帰式によって説明される部分。
y_i - \hat{y}_i
残差変動。平均からのばらつきのうち、回帰式によって説明されない部分。

当然、[全変動]＝[回帰変動]＋[残差変動]となることはすぐに確認できる。

また、これら3つの変動の2乗和をTSS(Total Sum of Squares)、ESS(Explained Sum of Squares、RSS(Residual Sum of Squares)として、以下のように定義する。

(8)   \begin{eqnarray*} TSS &=& \sum_{i=1}^n (y_i - \overline{Y})^2 \\ ESS &=& \sum_{i=1}^n (\hat{y}_i - \overline{Y})^2 \\ RSS &=& \sum_{i=1}^n (y_i - \hat{y}_i)^2 \end{eqnarray*}

このとき、それぞれについて以下のように表せる。

(9)   \begin{eqnarray*} TSS &=& n {\sigma _Y}^2 \\ ESS &=& n \frac{{\sigma_{XY}}^2}{{\sigma_x}^2} \\ RSS &=& n {\sigma_Y}^2 - n \frac{{\sigma_{XY}}^2}{{\sigma_X}^2} \end{eqnarray*}

これらから、全変動・回帰変動・残差変動の2乗和についても、以下の関係が成り立つことがわかる。

(10)   \begin{equation*} TSS = ESS + RSS \end{equation*}

補足～ESS

式(9) のESSの確認。以下のように、線形変換された分散から導ける。

(11)   \begin{eqnarray*} ESS &=& \sum_{i=1}^n \left( a x_i +b - (a \overline{X} + b) \right)^2 \\ &=& n {\rm Var}(a X + b) \\ &=& na^2{\rm Var}(X) \\ &=& n a^2 {\sigma_x}^2 \\ &=& n \left( \frac{\sigma_{XY}}{{\sigma_x}^2} \right)^2 {\sigma_x}^2 \\ &=& n \frac{{\sigma_{XY}}^2}{{\sigma_x}^2} \end{eqnarray*}

あるいは、以下のように地道に計算しても確認できる。

(12)   \begin{eqnarray*} \sum_{i=1}^n (\hat{y}_i - \overline{Y}) &=& \sum_{i=1}^n ({\hat{y}_i}^2 - 2 \hat{y}_i \overlne{Y} + \overline{Y}^2) \\ &=& \sum_{i=1}^n (\hat{y}_i^2 - 2(a x_i + b) \overline{Y} +\overline{Y}^2) \\ &=& \sum_{i=1}^n (\hat{y}_i^2 - 2a x_i \overline{Y} - 2b\overline{Y} + \overline{Y}^2) \\ &=& \sum_{i=1}^n \hat{y}_i^2 - 2an \overline{X} \; \overline{Y} -2nb \overline{Y} + n \overline{Y}\end{eqnarray*} \\

ここで

(13)   \begin{eqnarray*} \sum_{i=1}^n \hat{y}_i^2 &=& \sum_{i=1}^n (ax_i + b)^2 \\ &=&\sum_{i=1}^n (a^2 x_i^2 + 2ab x_i + b^2) \\ &=& a^2 \sum_{i=1}^n x_i^2 + 2abn \overline{X} + nb^2 \end{eqnarray*}

これを代入して、

(14)   \begin{eqnarray*} ESS &=& a^2 \sum_{i=1}^n x_i^2 + 2nab \overline{X} + nb^2 - 2na \overline{X} \; \overline{Y} - 2nb \overline{Y} + n \onerline{Y} \end{eqnarray*}

さらに\sum_{i=1}^n x_i^2 = \n \sigma_X^2 + n \overline{X}^2を考慮して、

(15)   \begin{eqnarray*} ESS &=& n a^2 {\sigma_X}^2 + n a^2 \overline{X}^2 + 2nab \overline{X} + nb^2 -2na\overline{X}\;\overline{Y} - 2nb\overline{Y} + n\overline{Y}^2 \\ &=& n a^2 {\sigma_X}^2 + n(a \overline{X} - \overline{Y})^2 + 2nab\overline{X} + nb^2 - 2nb\overline{Y} \\ &=& n a^2 {\sigma_X}^2 + nb^2 + 2nab\overline{X} + nb^2 - 2nb\overline{Y} \\ &=& n a^2 {\sigma_X}^2 + 2nb(b + a\overline{X} - \overline{Y}) \\ &=& n a^2 {\sigma_X}^2 \\ &=& n \left( \frac{\sigma_{XY}}{{\sigma_X}^2} \right)^2 {\sigma_X}^2 \\ &=& n \frac{{\sigma_{XY}}^2}{{\sigma_X}^2} \end{eqnarray*}

補足～RSS

式(9) のRSSの確認。

(16)   \begin{eqnarray*} RSS &=& \sum_{i=1}^n (y_i - \hat{y}_i)^2 = \sum_{i=1}^n (y_i - a x_i - b)^2 \\ &=& \sum_{i=1}^n (y_i - a x_i - \overline{Y} + a \overline{X})^2 \\ &=& \sum_{i=1}^n \left(y _i - \overline{Y} - a (x_i - \overline{X}) \right) ^2 \\ &=& \sum_{i=1}^n \left(y _i - \overline{Y} - \frac{\sigma_{XY}}{{\sigma_X}^2} (x_i - \overline{X}) \right) ^2 \\ &=& n {\sigma_Y}^2 - 2n \frac{\sigma_{XY}}{{\sigma_X}^2} \sigma_{XY} + n \frac{{\sigma_{XY}}^2}{{\sigma_X}^4} {\sigma_X}^2 \\ &=& n {\sigma_Y}^2 - 2n \frac{{\sigma_{XY}}^2}{{\sigma_X}^2} + n \frac{{\sigma_{XY}}^2}{{\sigma_X}^2} \\ &=& n {\sigma_Y}^2 - n \frac{{\sigma_{XY}}^2}{{\sigma_X}^2} \end{eqnarray*}

決定係数

決定係数は、全変動のうち回帰変動でどれだけ説明できるか（残差変動が少ないか）を表したもので、通常R^2で表される。

(17)   \begin{align*} R^2 = \frac{ESS}{TSS} = \frac{TSS - RSS}{TSS} = 1 - \frac{RSS}{TSS} \end{align*}

この値は、残差変動が小さいほど1に近づき、RSS = TSSすなわち全変動が回帰変動で全く説明できないときに0となる。