Advanced Biomedical Engineering

近藤の修論をまとめた論文がAdvanced Biomedical Engineeringにacceptされました.
T. Uchiyama, G. Kondo: Relationships among electromyogram, displacement and velocity of the center of pressure, and muscle stiffness of the medial gastrocnemius muscle during quiet standing

第59回日本生体医工学会大会

小倉の卒業研究をまとめて第59回日本生体医工学会大会(5月25日から27日)で発表しました.岡山で開催の予定でした.オンラインで開催されました.
T. Uchiyama and Y. Ogura: Accurate natural frequency estimation method for myotonometer using a system identification method

もう1件,連名で発表しました.
赤澤堅造,奥野竜平,内山孝憲:骨格筋の新しい数学モデル−フィラメント滑走SLモードと伸展による力増強STモード

IEEE EMBC 2019

ベルリンで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
直線は固有ベクトルの方向を表す
第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]
固有ベクトルの符号が反転していると左右・上下に反転がある
主成分分析のライブラリで計算したもの

第1主成分は,固有値が最も大きい固有ベクトルの方向です.固有値が最も大きい方向は,データの散らばりが最も大きくなる方向です.これは次のように導くことができます.


ベクトル\mathbf x_iの平均は\bar{\mathbf x}=\displaystyle\frac{1}{N}\sum_{i=1}^{N}\mathbf x_iです.
今,データのベクトル\mathbf x_iをデータの散らばりが最も大きくなる方向のベクトル\mathbf u_1に射影することを考えます.このとき,\mathbf u_1を単位ベクトルにとります.\mathbf x_i\mathbf u_1に射影すると\mathbf x_i^T\mathbf u_1になります.射影の分散は
\sigma^2 = \displaystyle\frac{1}{N}\sum_{i=1}^{N}(\mathbf x_i^T\mathbf u_1 - \bar{\mathbf x}^T\mathbf u_1)^2 = \mathbf u_1^T S \mathbf u_1です.
ここで
S= \displaystyle\frac{1}{N}(\mathbf x_i - \bar{\mathbf x})(\mathbf x_i - \bar{\mathbf x})^Tです.
\mathbf u_1^T S \mathbf u_1\mathbf u_1^T\mathbf u_1=1の条件の下で最大化します.これはラグランジェの未定乗数法を用いて
\mathbf u_1^T S \mathbf u_1 - \lambda(\mathbf u_1^T\mathbf u_1 - 1)
を最大化して解きます.\mathbf u_1で微分すると
S\mathbf u_1 - \lambda \mathbf u_1=0
になります.つまりS\mathbf u_1 = \lambda \mathbf u_1ですから,\mathbf u_1は行列Sの固有ベクトルです.
この式に左から\mathbf u_1^Tをかけると\mathbf u_1^T S \mathbf u_1 = \lambdaとなり,データ\mathbf x_iの射影の分散は固有値です.また,それを最大にするものは最大の固有値になります.