多変量解析入門(isbn:4884122801 篠原出版新社 足立堅一)の「13.3.2 変動の分解」の
ST = SS + CF
の式が成立することをsympyを使って実際にやってみることにします。
今回はisympyを使います。まずはシンボルの準備
In [1]: x1 = Symbol('x1')
In [2]: x2 = Symbol('x2')
In [3]: x3 = Symbol('x3')
次にSTを設定
In [4]: ST = x1**2 + x2**2 + x3**2 In [5]: ST Out[5]: 2 2 2 x₁ + x₂ + x₃
SSを設定する前にmを設定
In [6]: m = (x1 + x2 + x3)/3 In [7]: m Out[7]: x₁ x₂ x₃ ── + ── + ── 3 3 3
SSをmを使って設定
In [8]: SS = (x1-m)**2 + (x2-m)**2 + (x3-m)**2
In [9]: SS
Out[9]:
2 2 2
⎛ x₁ x₂ 2*x₃⎞ ⎛ x₂ x₃ 2*x₁⎞ ⎛ x₁ x₃ 2*x₂ ⎞
⎜- ── - ── + ───⎟ + ⎜- ── - ── + ───⎟ + ⎜- ── - ── + ───⎟
⎝ 3 3 3 ⎠ ⎝ 3 3 3 ⎠ ⎝ 3 3 3 ⎠
まず素直にCFを(CF0として)求める
In [10]: CF0 = ST - SS In [11]: CF0 Out[11]: （長いので省略）
expand()をつかって単純化(CF1とする)。
In [12]: CF1 = CF0.expand()
In [13]: CF1
Out[13]:
2 2 2
2*x₁*x₂ 2*x₁*x₃ 2*x₂*x₃ x₁ x₂ x₃
─────── + ─────── + ─────── + ─── + ─── + ───
3 3 3 3 3 3
ここで、”k = x1 + x2 + x3”という式をaとして設定
In [15]: k = Symbol('k')
In [16]: a = Eq(k, x1 + x2 + x3)
In [17]: a
Out[17]: k = x₁ + x₂ + x₃
aの式をx1を求める形式に変更し、それをCF1に当てはめてCF2とする。”[0]”はsolveがlistを返すため。
In [18]: solve(a, x1) Out[18]: [k - x₂ - x₃] In [19]: CF2 = CF1.subs(x1, solve(a, x1)[0]) In [20]: CF2 Out[20]: （長いので省略）
CF2をexpand()してCF3にする
In [21]: CF3 = CF2.expand() In [22]: CF3 Out[22]: 2 k ── 3
ラスト。aを使ってkを消す。
In [23]: CF3.subs(k, solve(a, k)[0])
Out[23]:
2
(x₁ + x₂ + x₃)
───────────────
3
できました。
さらっと書いているけど、できるようになるまでかなり時間がかかりました。
（それにしても、小さい数字がコピペで入ったのはビックリ。これWindowsでみても綺麗にでるのかな？）
