![](https://www.bi.appi.keio.ac.jp/wordpress/wp-content/uploads/2019/11/IMG_3127-300x225.jpg)
堀研究室と合同で四万温泉へ行きました.写真は四万湖です.
物理情報数学B,バイオシステムおよび実験レポートのグラフについて,Pythonのプログラムの例へのリンクのページを作成しました.1, 2, 3年生の皆さんへのページにリンクがあります.
ベルリンで7月23日から27日まで開催された41st Annual International Conference of the IEEE Engineering in Medicine and Biology Societyで
T. Uchiyama, E. Hamada, Relationship between Medial Gastrocnemius Muscle Stiffness and the Angle between the Rearfoot and Floor
を発表しました.
足関節角度と腓腹筋および足関節のスティフネスの関係(かかとの高さがある靴を履いたときのように,足の母指球より後側を斜面に置いて立っている状態です)を調べました.足関節角度が増加するにつれて腓腹筋のスティフネスは増加しました.
物理情報数学Bで扱った行列AをA=UΣVTに特異値分解します.
# import numpy as np import matplotlib.pyplot as plt A = np.array([[2, 2], [-1, 2]]) U, S, V = np.linalg.svd(A) # V^TがVに入ることに注意する print("特異値分解") print("U") print(U) print("S") print(S) print("V") print(V) print("USV^T") S2 = np.eye(2)*S # plt.figure(figsize=(5, 5)) plt.xlim([-3, 3]) plt.ylim([-3, 3]) plt.xlabel('x1') plt.ylabel('x2') print(np.dot(np.dot(U, S2), V)) # theta = np.arange(0, 361, 3) uc = np.array([np.cos(theta/180*np.pi), np.sin(theta/180*np.pi)]) # print("Vの列ベクトルと単位円") plt.plot(uc[0,:], uc[1,:], color='black', alpha=0.3) plt.quiver(0, 0, V[0, 0], V[0, 1], angles="xy", scale_units='xy', scale=1, color='red') plt.quiver(0, 0, V[1, 0], V[1, 1], angles="xy", scale_units='xy', scale=1, color='blue') plt.show() x1 = np.dot(V, V.T) c1 = np.dot(V, uc) print("Vの列をV^Tで写像する V^TV 回転と鏡映") print(x1) plt.figure(figsize=(5, 5)) plt.xlim([-3, 3]) plt.ylim([-3, 3]) plt.xlabel('x1') plt.ylabel('x2') plt.plot(c1[0,:], c1[1,:], color='black', alpha=0.3) plt.quiver(0, 0, x1[0, 0], x1[1, 0], angles="xy", scale_units='xy', scale=1, color='red') plt.quiver(0, 0, x1[0, 1], x1[1, 1], angles="xy", scale_units='xy', scale=1, color='blue') plt.show() print("それをシグマで写像する ΣV^TV 軸方向に拡大・縮小(特異値倍)する") x2 = np.dot(S2, x1 ) c2 = np.dot(S2, c1) print(x2) plt.figure(figsize=(5, 5)) plt.xlim([-3, 3]) plt.ylim([-3, 3]) plt.xlabel('x1') plt.ylabel('x2') plt.plot(c2[0,:], c2[1,:], color='black', alpha=0.3) plt.quiver(0, 0, x2[0, 0], x2[1, 0], angles="xy", scale_units='xy', scale=1, color='red') plt.quiver(0, 0, x2[0, 1], x2[1, 1], angles="xy", scale_units='xy', scale=1, color='blue') plt.show() print("最後にUで写像する UΣV^TV 回転と鏡映") print("u1とu2のσ1倍とσ2倍にそれぞれ重なる") x3 = np.dot(U, x2) c3 = np.dot(U, c2) print(x3) plt.figure(figsize=(5, 5)) plt.xlim([-3, 3]) plt.ylim([-3, 3]) plt.xlabel('x1') plt.ylabel('x2') plt.plot(c3[0,:], c3[1,:], color='black', alpha=0.3) plt.quiver(0, 0, x3[0, 0], x3[1, 0], angles="xy", scale_units='xy', scale=1, color='red') plt.quiver(0, 0, x3[0, 1], x3[1, 1], angles="xy", scale_units='xy', scale=1, color='blue') plt.quiver(0, 0, S[0]*U[0, 0], S[0]*U[1, 0], angles="xy", scale_units='xy', scale=1, color='black', alpha=0.3) plt.quiver(0, 0, S[1]*U[0, 1], S[1]*U[1, 1], angles="xy", scale_units='xy', scale=1, color='black', alpha=0.3) plt.show()
特異値分解 U [[ 0.89442719 0.4472136 ] [ 0.4472136 -0.89442719]] S [3. 2.] V [[ 0.4472136 0.89442719] [ 0.89442719 -0.4472136 ]] USV^T [[ 2. 2.] [-1. 2.]] Vの列ベクトルと単位円
Vの列をV^Tで写像する V^TV 回転と鏡映 [[ 1.00000000e+00 -2.43158597e-17] [-2.43158597e-17 1.00000000e+00]]
それをシグマで写像する ΣV^TV 軸方向に拡大・縮小(特異値倍)する [[ 3.00000000e+00 -7.29475792e-17] [-4.86317195e-17 2.00000000e+00]]
最後にUで写像する UΣV^TV 回転と鏡映 u1とu2のσ1倍とσ2倍にそれぞれ重なる [[ 2.68328157 0.89442719] [ 1.34164079 -1.78885438]]
物理情報数学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
主成分分析のライブラリを使う 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] 固有ベクトルの符号が反転していると左右・上下に反転がある
第1主成分は,固有値が最も大きい固有ベクトルの方向です.固有値が最も大きい方向は,データの散らばりが最も大きくなる方向です.これは次のように導くことができます.
ベクトルの平均は
です.
今,データのベクトルをデータの散らばりが最も大きくなる方向のベクトル
に射影することを考えます.このとき,
を単位ベクトルにとります.
を
に射影すると
になります.射影の分散は
です.
ここでです.
を
の条件の下で最大化します.これはラグランジェの未定乗数法を用いて
を最大化して解きます.で微分すると
になります.つまりですから,
は行列
の固有ベクトルです.
この式に左からをかけると
となり,データ
の射影の分散は固有値です.また,それを最大にするものは最大の固有値になります.
物理情報数学Bで対称行列の応用例として取り上げている傾いた楕円の式の説明です.
の式を行列で表すと
になる.
行列をと置くと,
の固有ベクトル(正規化されたもの)を列に持つ行列
を用いて
と表せる.行列
は固有値として9と1を持つ.固有値が9のとき固有ベクトルは
である.固有値が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が重いので結果だけ載せます.
で固有値の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]]
で固有値の1つが正で1つが負の場合
固有値 [ 1. -1.] 固有ベクトル [[ 0.70710678 0.70710678] [-0.70710678 0.70710678]]
3年生秋学期後半の選択科目,バイオシステムで扱う内容の一部をPythonで実装したものです.keio IDが必要です.
線形代数の計算をPythonのプログラムで表しました.keio IDが必要です.
Pythonで簡単なグラフを描くためのサンプルです.Google Colaboratoryを使います.keio IDが必要です.