主成分分析

物理情報数学Bで説明していている対称行列が出てくる例です.固有ベクトルへの射影も関係します.

#
# 主成分分析
#
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

N = 100
d = 16
x = 50 + np.random.randn(N) * d
y = 0.8 * x + np.random.randn(N) * d * 0.5 + 10
plt.figure(figsize=(5, 5))
plt.xlim([0, 100])
plt.ylim([0, 100])
plt.xlabel('x')
plt.ylabel('y')
plt.scatter(x, y, color='blue', alpha=0.5)

features = np.vstack([x, y]).T
features2 = features.copy()
ave = np.average(features,0)
print("平均")
print(ave)
features = features - ave
cv = np.dot(features.T, features) / (N - 1)
print("共分散行列")
print(cv)
W, v = np.linalg.eig(cv)
print("共分散行列の固有ベクトルを列にもつ行列")
print("{0}".format(v))
print("固有ベクトル")
print("{0}".format(v[:,0].reshape(2,1)))
print("固有ベクトル")
print("{0}".format(v[:,1].reshape(2,1)))
print("共分散行列の固有値:{0}".format(W))
# 寄与率は固有値/(固有値の総和)
p = np.trace(cv)
print("寄与率1: {0}, 2: {1}".format(W[0]/p, W[1]/p))
# 固有値の大きい順にソート
ind = np.argsort(W)[::-1]
sortedW = np.sort(W)[::-1]
sortedv = v[:,ind]
# 固有ベクトル方向の直線
plt.plot(np.array([0, 100]), np.array([sortedv[1,0]/sortedv[0,0]*(0 - ave[0])+ave[1], sortedv[1,0]/sortedv[0,0]*(100 - ave[0])+ave[1]]),  color='blue', alpha=0.5)
plt.plot(np.array([0, 100]), np.array([sortedv[1,1]/sortedv[0,1]*(0 - ave[0])+ave[1], sortedv[1,1]/sortedv[0,1]*(100 - ave[0])+ave[1]]),  color='blue', alpha=0.5)
#
t1 = np.zeros([2, N])
# 固有ベクトル方向に射影する
#p1 = v[:,0].reshape(2,1)
#p2 = v[:,1].reshape(2,1)
#for i in range(N):
#  t1[0, i] = np.dot(p1.T, (features.T[:,i]).reshape(2,1)) # 射影するベクトルは単位ベクトルなので2乗ノルムの2乗で割る必要なし
#  t1[1, i] = np.dot(p2.T, (features.T[:,i]).reshape(2,1))
t1 = np.dot(sortedv.T, features.T)
#
plt.figure(figsize=(5, 5))
plt.xlim([-70, 70])
plt.ylim([-70, 70])
plt.xlabel('p1')
plt.ylabel('p2')
plt.scatter(t1[0,:], t1[1,:], color='red', alpha=0.5)
# 3σの楕円
theta = np.arange(0, 361, 3)
X = np.cos(theta/180*np.pi) * np.sqrt(sortedW[0]) * 3
Y = np.sin(theta/180*np.pi) * np.sqrt(sortedW[1]) * 3
plt.plot(X, Y, color='red', alpha=0.8)
plt.show()
#
#
print("主成分分析のライブラリを使う")
pca = PCA(n_components=2)
pca.fit(features2)
print("components")
print("{0}".format(pca.components_))
print("mean={0}".format(pca.mean_))
print("covariance")
print("{0}".format(pca.get_covariance()))
t2 = pca.fit_transform(features)
print("contribution")
print(pca.explained_variance_ratio_)
plt.figure(figsize=(5, 5))
plt.xlim([-70, 70])
plt.ylim([-70, 70])
plt.xlabel('p1')
plt.ylabel('p2')
plt.scatter(t2[:,0], t2[:,1], color='green', alpha=0.5)
print("固有ベクトルの符号が反転していると左右・上下に反転がある")
plt.show()
平均
[52.53235793 51.86342693]
共分散行列
[[244.04297958 193.42743647]
 [193.42743647 216.02097448]]
共分散行列の固有ベクトルを列にもつ行列
[[ 0.73220426 -0.6810851 ]
 [ 0.6810851   0.73220426]]
固有ベクトル
[[0.73220426]
 [0.6810851 ]]
固有ベクトル
[[-0.6810851 ]
 [ 0.73220426]]
共分散行列の固有値:[423.96619622  36.09775785]
寄与率1: 0.9215375220555805, 2: 0.07846247794441952
直線は固有ベクトルの方向を表す
第1主成分と第2主成分を軸にとったもの.3σの長半径と短半径を持つ楕円を重ねて描いた
主成分分析のライブラリを使う
components
[[ 0.73220426  0.6810851 ]
 [-0.6810851   0.73220426]]
mean=[52.53235793 51.86342693]
covariance
[[244.04297958 193.42743647]
 [193.42743647 216.02097448]]
contribution
[0.92153752 0.07846248]
固有ベクトルの符号が反転していると左右・上下に反転がある
主成分分析のライブラリで計算したもの

傾いた楕円の式

物理情報数学Bで対称行列の応用例として取り上げている傾いた楕円の式の説明です.

5x^2 + 8xy + 5y^2 = 1
の式を行列で表すと
\left[\begin{array}{cc}x & y\end{array}\right]\left[\begin{array}{cc}5 & 4\\ 4 & 5\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]=1
になる.

行列をAと置くと,Aの固有ベクトル(正規化されたもの)を列に持つ行列Qを用いてA=Q\Lambda Q^Tと表せる.行列Aは固有値として9と1を持つ.固有値が9のとき固有ベクトルは[1\ 1]^Tである.固有値が1のとき固有ベクトルは[1\ -1]^Tである.固有ベクトルの大きさは\sqrt{2}であるので
\left[\begin{array}{cc}x & y\end{array}\right]\displaystyle\frac{1}{\sqrt{2}}\left[\begin{array}{rr}1 & 1\\ 1 & -1\end{array}\right]\left[\begin{array}{cc}9 & 0\\ 0 & 1\end{array}\right]\frac{1}{\sqrt{2}}\left[\begin{array}{rr}1 & 1\\ 1 & -1\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]=1
となる.

この式は
9\left(\displaystyle\frac{x+y}{\sqrt{2}}\right)^2 + 1\left(\displaystyle\frac{x-y}{\sqrt{2}}\right)^2=1
と表せる.括弧内をXYと置くと,
9X^2+Y^2=1
の楕円の式になる.楕円の短半径は1/3で長半径は1になる.これらは,固有値の逆数の平方根である.

import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg as LA

#
# 5x^2 + 8xy + 5y^2 = 1
# 9X^2 + y^2 = 1
# X = (x + y)/sqrt(2)
# Y = (x - y) /sqrt(2)
#
theta = np.arange(0, 361, 3)
X = np.cos(theta/180*np.pi)/3
Y = np.sin(theta/180*np.pi)
x = (X+Y)/np.sqrt(2)
y = (X-Y)/np.sqrt(2)
plt.figure(figsize=(5, 5))
plt.xlim([-1.2, 1.2])
plt.ylim([-1.2, 1.2])
plt.xlabel('x')
plt.ylabel('y')
plt.plot(x, y, color='blue', alpha=0.8)
#
A = np.array([[5, 4], [4, 5]])
w,v = LA.eig(A)
print('固有値')
print(w)
print('固有ベクトル')
print(v)
# 固有値の逆数の平方根が長半径と短半径
# 固有ベクトルが軸方向
a = v[:,0].reshape(2, 1) / np.sqrt(w[0])
b = v[:,1].reshape(2 ,1) / np.sqrt(w[1])
plt.plot(np.array([0, a[0]]), np.array([0, a[1]]), color='red')
plt.plot(np.array([0, b[0]]), np.array([0, b[1]]), color='green')
plt.show()
固有値
[9. 1.]
固有ベクトル
[[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]

固有値・固有ベクトルと行列微分方程式

Google Colaboratoryが重いので結果だけ載せます.

\displaystyle\frac{dx_1(t)}{dt} = x_2(t)\\ \displaystyle\frac{dx_2(t)}{dt}=x_1(t)
で固有値の1つが正で1つが負の場合

import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg as LA
import random

A = np.array([[0, 1], [1, 0]])
w, v = LA.eig(A)
print('固有値')
print(w)
print('固有ベクトル')
print(v)
x10 = 4
x20 = 2
t = np.arange(0.0, 2.0, 0.1)
x10 = [random.uniform(-10, 10) for i in range(300)]
x20 = [random.uniform(-10, 10) for i in range(300)]

N = 22
x = np.linspace(-10,  10,  N)
y = np.linspace(-10,  10,  N)
(xm, ym) = np.meshgrid(x, y)
dxdt = ym
dydt = xm
norm = (np.sqrt(dxdt**2 + dydt**2))
dxdt2 = dxdt/norm
dydt2 = dydt/norm
plt.figure(figsize=(5, 5))
plt.quiver(xm, ym, dxdt2, dydt2, angles='xy', scale_units='xy', scale=1)
plt.xlim([-10, 10])
plt.ylim([-10, 10])
plt.xlabel('x1')
plt.ylabel('x2')
plt.show()

# 解析解
plt.figure(figsize=(5, 5))
for i in range(300):
  c = LA.solve(v, np.array([x10[i], x20[i]]))
  xxa = c[0] * np.exp(w[0]*t) * v[:,0].reshape(2, 1) +c[1] * np.exp(w[1]*t) * v[:,1].reshape(2, 1)
  plt.plot(xxa[0,:], xxa[1,:], color='blue', alpha=0.2)
plt.xlim([-10, 10])
plt.ylim([-10, 10])
plt.xlabel('x1')
plt.ylabel('x2')
plt.show()
固有値
[ 1. -1.]
固有ベクトル
[[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]

\displaystyle\frac{dx_1(t)}{dt} = -x_2(t)\\ \displaystyle\frac{dx_2(t)}{dt}=-x_1(t)
で固有値の1つが正で1つが負の場合

固有値
[ 1. -1.]
固有ベクトル
[[ 0.70710678  0.70710678]
 [-0.70710678  0.70710678]]

餌と捕食者の関係

物理情報数学Bで,行列の累乗の応用例として,餌と捕食者の関係を表す行列差分方程式を説明しています.

「第k年のイワシの数をx_1(k),クジラの数をx_2(k)とする.クジラがいないときの1年後のイワシの数を\alpha x_1(k)とする.つまり,第k年のイワシの数に比例して変化する.クジラに1年間に食べられるイワシの数を\beta x_2(k)とする.つまり,クジラの数に比例するイワシが食べられる.イワシがいないときの1年後のクジラの数を\gamma x_2(k)とする.また,イワシを食べて1年間に増加するクジラの数を\delta x_1(k)とする.」

上記の関係を行列差分方程式で表すと

   \left[\begin{array}{c}x_1(k+1) \\x_2(k+2)\end{array}\right]   =   \left[\begin{array}{rr}\alpha & -\beta \\ \delta & \gamma\end{array}\right]   \left[\begin{array}{c}x_1(k) \\x_2(k)\end{array}\right]

になります.一般解は

\boldsymbol{x}(k) = \left[\begin{array}{rr}\alpha & -\beta \\ \delta & \gamma\end{array}\right]^k\boldsymbol{x}(0)

になります.

この行列差分方程式を講義では手計算で解けるように,\alpha=1.5\beta=0.25\gamma=1.1\delta=0.07として解きます.しかし,この行列の固有値は1以上の実数が2つですので,餌と捕食者の数の変化に期待される振動的な振る舞いにはなりません.初期値のクジラの数が多いと,下図のようにイワシの数は減少に転じます.

Graph0
餌と捕食者の数(x1(0)=100,x2(0)=150)

餌と捕食者の関係のモデルとしては,ロトカ・ボルテラのモデルが知られています.このモデルでは,

   \displaystyle\frac{dx_1(t)}{dt} = \alpha x_1(t) - \beta x_1(t)x_2(t) \\   \displaystyle\frac{dx_2(t)}{dt} = \delta x_1(t)x_2(t) - \gamma x_2(t)

で表されます.非線形微分方程式ですから,x_1(t)x_2(t)を解析的に求めることはできません.適当な値を代入して数値解を求めると,x_1(t)x_2(t)が振動的な振る舞いをすることがわかります.

また,平衡点はdx_1(t)/dt=0dx_2(t)/dt=0から,原点とx_1=\gamma/\deltax_2=\alpha/\betaにあります.平衡点近傍では,x_1(t)x_2(t)の項が小さいので無視できると仮定すると,(0, 0)の平衡点近傍では,

   \displaystyle\frac{d}{dt}\left[\begin{array}{c}x_1(t) \\x_2(t)\end{array}\right]   =   \left[\begin{array}{rr}\alpha & 0 \\   0 & -\gamma\end{array}\right]   \left[\begin{array}{c}x_1(t) \\x_2(t)\end{array}\right]

となり,行列の固有値は実数で1つは正でもう1つは負です.したがって,(0, 0)は鞍点です.

また,もう1つの平衡点(\gamma/\delta, \alpha/\beta)の近傍では,x_1(t)x_1+\gamma/\deltax_2(t)x_2+\alpha/\betaとおいて,x_1(t)x_2(t)の項が小さいので無視できると仮定すると

   \displaystyle\frac{d}{dt}\left[\begin{array}{c}x_1(t) \\x_2(t)\end{array}\right]   =   \left[\begin{array}{rr}0 & -\frac{\beta\gamma}{\delta} \\   \frac{\alpha\delta}{\beta} & 0\end{array}\right]   \left[\begin{array}{c}x_1(t) \\x_2(t)\end{array}\right]

と表されます.したがって,固有値が純虚数となり,振動する解になります.

\alpha=10\beta=0.1\gamma=10\delta=0.1とし,x_1(0)=10x_2(0)=100 で数値解を求めると下図のようになります.上のパネルはx_1(t)x_2(t)の時間発展です. 下のパネルはx_1(t)x_2(t)の関係です.

x1(t)とx2(t)の時間発展

lv2
x1(t)とx2(t)の関係(t=0で(10, 100)とx1=0のヌルクラインに近く,x2が急速に減少するが,原点が鞍点であるのである程度減少したところでx1が急速に増加する.)

初期値を平衡点((100, 100))に近い点(99, 99)とすると下図のように円に近い軌道を描きます.

lv3
初期値を(99, 99)と平衡点(100, 100)に近い値にしたので円に近い軌道を描く

MATLABのコードです.

  • lv.m
tspan = [0:0.001:5];
[t y] = ode45(@LotkaVolterra,tspan, [99 99]);
figure(1);
plot(t, y(:,1), 'b', t, y(:,2), 'r');
figure(2);
plot(y(:,1), y(:,2));
  • LotkaVolterra.m
function dx = LotkaVolterra(t, x)
alpha = 10;
beta = 0.1;
gamma = 10;
delta = 0.1;
dx = zeros(2,1);
dx(1) = alpha * x(1) - beta * x(1) * x(2);
dx(2) = delta * x(1) * x(2) - gamma * x(2);
end

自然科学実験(物理学)オシロスコープ

RC回路の周波数特性を計測する実験のレポートで,理論値と実験値のずれに関する考察でオシロスコープの内部抵抗を原因の一つに上げているものがしばしばあります.結論から言えば,実験の回路ではオシロスコープの内部抵抗の影響は小さいです.

オシロスコープは電圧計です.理想的な電圧計では,電流が全く電圧計に流れないのですが,現実的には流れ込むため,理想的な電圧計とそれに並列な有限の抵抗で,現実の電圧計の等価回路を表します.オシロスコープの内部インピーダンスは1 MΩ程度です(プローブを使わないでリード線で配線しているので,10:1プローブを付けているときのように10 MΩ程度にはなりません).それに対して実験に使っている回路の抵抗は10 kΩです.1/100ですから,影響は小さいです.むしろ抵抗やコンデンサの容量の誤差の方が大きいでしょう.

_MG_3600

ハイパスフィルタ

自然科学実験(物理学)のテーマ,オシロスコープの応用課題で特性を調べるローパスフィルタの抵抗とコンデンサを入れ替えて構成したハイパスフィルタです.

rc-highpass

 

v_i(t) = \displaystyle\frac{1}{C}\int i(t) dt + v_o(t)v_o(t) = R i(t)の関係がありますから,

v_i(t) = \displaystyle\frac{1}{RC}\int v_o(t) dt + v_o(t)

になります.

ラプラス変換で解く

ラプラス変換して,

V_i(s) = \displaystyle\frac{V_o(s)}{RCs} + V_o(s)

を得ます.伝達関数は

\displaystyle\frac{V_o(s)}{V_i(s)} = \frac{1}{\frac{1}{RCs} + 1} = \frac{s}{s+\frac{1}{RC}}

になります.伝達関数の極は,抵抗とコンデンサを入れ替えて構成するローパスフィルタと同じになります.形式的にs=j\omegaとすると,ゲイン特性は

\left|\displaystyle\frac{V_o(s)}{V_i(s)}\right| = \displaystyle\frac{\omega}{\sqrt{\omega^2 + \left(\frac{1}{RC}\right)^2}}

になります.

位相特性は

90^\circ - \tan^{-1}\omega RC

になります.

1年生の数学の範囲で解く

同次微分方程式の解,つまり基本解は時間とともに零になるので,v_i(t)=\sin\omega tのときの特解を求めればよいのは,ローパスフィルタの場合と同様です.

積分を含む式

v_i(t) = \displaystyle\frac{1}{RC}\int v_o(t) dt + v_o(t)

なので,tで微分して,

\displaystyle\frac{dv_i(t)}{dt}=\frac{1}{RC}v_o(t) + \frac{dv_o(t)}{dt}

を得ます.\displaystyle\frac{dv_i(t)}{dt}=\omega\cos\omega tを代入すると,

\omega\cos\omega t=\displaystyle\frac{1}{RC}v_o(t) + \frac{dv_o(t)}{dt}

になります.v_o(t)は,v_i(t)と同じ角周波数の正弦波ですが,振幅と位相が異なるものになるので,解をa\sin(\omega t + \phi)とおいて,a\phi を求めます.

\omega\cos\omega t = \displaystyle\frac{1}{RC}a\sin(\omega t + \phi) + a\omega\cos(\omega t + \phi)

になるので三角関数の合成を用いて,

\omega\cos\omega t = a\sqrt{\displaystyle\left(\frac{1}{RC}\right)^2 + \omega^2}\sin(\omega t + \phi + \psi)

\psi=\tan^{-1}RC\omega

となり,両辺を\cosにして

\displaystyle\omega\cos\omega t = a\sqrt{\displaystyle\left(\frac{1}{RC}\right)^2 + \omega^2}\cos\left(-\frac{\pi}{2}+\omega t + \phi + \psi\right)

となります.したがって

a=\displaystyle\frac{\omega}{\sqrt{\left(\frac{1}{RC}\right)^2 + \omega^2}}

を得ます.

位相は\phi=\displaystyle\frac{\pi}{2}-\psiなので,

90^\circ - \tan^{-1}RC\omega

になります.

自然科学実験(物理学) オシロスコープ

1年生の自然科学実験のテーマの一つにオシロスコープがあります.応用課題として,抵抗とコンデンサで作るローパスフィルタの周波数特性を計測する実験があります.1年生では電気回路の講義がなく,交流回路の解析はまだ習っていません.1年生の数学では非同次微分方程式の解き方を学ぶので,タイミングによっては数学的な説明ができます.それすらできないときには,コンデンサの性質を定性的に説明することになります.

rc-lowpass

図の回路において,流れる電流をi(t)とおくと,抵抗Rの部分の電圧はRi(t)です.コンデンサCの電圧はv_o(t)です.また,i(t)=\displaystyle\frac{dv_o(t)}{dt}です(この式は,Q=CVQ=\displaystyle\int i(t)dtの関係で誘導しますが,高校の範囲では一定電流を\Delta t流すと電荷がQになることまでで,積分では習っていないと思います).

v_i(t)=Ri(t) + v_o(t)より,微分方程式v_i(t)=RC\displaystyle\frac{dv_o(t)}{dt} + v_o(t)を得ます.

1年生の数学の範囲で解く

1年生の数学で,非同次微分方程式の解き方を学んでいれば,v_i(t)=0のときの基本解は時間とともに減衰して零になるので,v_i(t)=\sin \omega tのときの特解を求めればよいことが分ります.

特解は,この回路は線形系ですから,入力の正弦波と同じ周波数で振動するすけれども,振幅と位相が異なるだけになります.そこで解をa\sin(\omega t + \phi)として,a\phiを求めることにします.

微分方程式に代入すると,

\sin\omega t=RCa\omega\cos(\omega t + \phi) + a\sin(\omega t + \phi)

になります.右辺に三角関数の合成を適用すると,

\sin\omega t = \sqrt{a^2 + (RCa\omega)^2}\sin(\omega t + \phi + \psi)

になります.両辺の振幅を比較して

a=\displaystyle\frac{1}{\sqrt{1+(RC\omega)^2}}

と,位相を比較して\phi = -\psiから

\phi = -\tan^{-1}RC\omega

を得ます.

交流回路の解析(フェーザ法)で解く

物理情報工学科では2年生春学期の電気回路同演習でフェーザ法を学びます.フェーザ法では,図の回路は

\dot{V_i} = R\dot{I} + \dot{V_o}

と表すことができ,

\dot{V_o} =\displaystyle\frac{1}{j\omega C}\dot{I}

を用いて\dot{I}を消去すると,

\dot{V_i} = (j\omega RC + 1)\dot{V_o}

を得ます.ゲイン特性は

\displaystyle\left|\frac{\dot{V_o}}{\dot{V_i}}\right| = \frac{1}{\sqrt{1 + (RC\omega)^2}}

になります.

位相特性は

-\tan^{-1}RC\omega

になります.

ラプラス変換で解く

物理情報工学科では2年生秋学期の物理情報数学Cでラプラス変換を学びます.ラプラス変換を学んでいれば,伝達関数に形式的にs=j\omegaとしてゲイン特性と位相特性を求めることができます.

微分方程式をラプラス変換すると

V_o(s) = RCsV_i(s) + V_i(s)

を得ます.伝達関数

\displaystyle\frac{V_o(s)}{V_i(s)} = \frac{1}{RCs + 1}

を得ます.s=j\omegaとすれば,上記と同じゲイン特性と位相特性になります.

補足

実験を早く終わった学生には,入力波形を矩形波にして積分されている様子を観察することや,抵抗とコンデンサを入れ替えてハイパスフィルタにして,定性的にですが,ハイパスフィルタの特性や入力波(三角波)が微分されている様子を観察してもらっています.