はじめに

　『スタンフォード　ベクトル・行列からはじめる最適化数学』の学習ノートです。
　「数式の行間埋め」や「Pythonを使っての再現」によって理解を目指します。本と一緒に読んでください。

　この記事は3.4節「角度」の内容です。
　3次元ベクトルのなす角と相関係数を比較するグラフを作成します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
2つのベクトルのなす角と相関係数の関係の可視化
参考書籍
おわりに

2つのベクトルのなす角と相関係数の関係の可視化

　2つの平均除去ベクトルのなす角(angle between two vectors)と相関係数(correlation coefficient)の関係をグラフで確認します。
　なす角については「【Python】3.4：2つのベクトルのなす角の可視化【『スタンフォード線形代数入門』のノート】 - からっぽのしょこ」、相関係数については「【Python】3.4：相関係数の可視化【『スタンフォード線形代数入門』のノート】 - からっぽのしょこ」を参照してください。
　「相関係数の可視化」では、2つのn次元ベクトルを散布図として、平面上で相関係数を可視化しました。この記事では、2つの3次元ベクトルを使って、なす角の大小と相関係数の大小などの関係から、ベクトルの類似度を可視化します。

　利用するライブラリを読み込みます。

# 利用ライブラリ
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

　3次元空間上のなす角の作図に使う配列の作成処理を関数として定義しておきます。

・コード(クリックで展開)

# なす角の作図用の処理を定義
def angle3d(a, b, alpha=1.0, beta=1.0):
    
    # xy面・xz面上のベクトルaのなす角を計算
    sgn_y = 1.0 if alpha*a[1] >= 0.0 else -1.0
    theta_a_xy = sgn_y * np.arccos(alpha*a[0] / np.linalg.norm(alpha * a[[0, 1]]))
    sgn_z = 1.0 if alpha*a[2] >= 0.0 else -1.0
    theta_a_xz = sgn_z * np.arccos(alpha*a[0] / np.linalg.norm(alpha * a[[0, 2]]))
    
    # xy面・xz面上のベクトルbのなす角を計算
    sgn_y = 1.0 if beta*b[1] >= 0.0 else -1.0
    theta_b_xy = sgn_y * np.arccos(beta*b[0] / np.linalg.norm(beta * b[[0, 1]]))
    sgn_z = 1.0 if beta*b[2] >= 0.0 else -1.0
    theta_b_xz = sgn_z * np.arccos(beta*b[0] / np.linalg.norm(beta * b[[0, 2]]))
    
    # 差が180°を超える場合は、負の角度を優角に変換
    if abs(theta_a_xy - theta_b_xy) > np.pi:
        theta_a_xy = theta_a_xy if alpha*a[1] >= 0.0 else theta_a_xy + 2.0*np.pi
        theta_b_xy = theta_b_xy if beta*b[1] >= 0.0 else theta_b_xy + 2.0*np.pi
    if abs(theta_a_xz - theta_b_xz) > np.pi:
        theta_a_xz = theta_a_xz if alpha*a[2] >= 0.0 else theta_a_xz + 2.0*np.pi
        theta_b_xz = theta_b_xz if beta*b[2] >= 0.0 else theta_b_xz + 2.0*np.pi
    
    # xy・xz面面上の角マーク用の値(ラジアン)を作成
    rad_xy_vals = np.linspace(start=theta_a_xy, stop=theta_b_xy, num=100)
    rad_xz_vals = np.linspace(start=theta_a_xz, stop=theta_b_xz, num=100)

    # xy面・xz面上の角マーク用の座標を計算
    r = 0.2
    angle_xy_x_vals = np.cos(rad_xy_vals) * r
    angle_xy_y_vals = np.sin(rad_xy_vals) * r
    angle_xz_x_vals = np.cos(rad_xz_vals) * r
    angle_xz_z_vals = np.sin(rad_xz_vals) * r
    angle_mark_xy = (angle_xy_x_vals, angle_xy_y_vals)
    angle_mark_xz = (angle_xz_x_vals, angle_xz_z_vals)
    
    # 角マーク用の角aの値(球面座標)を計算
    psi_a = np.arccos(alpha*a[2] / np.linalg.norm(alpha * a))
    sgn_y = 1.0 if alpha*a[1] >= 0.0 else -1.0
    phi_a = sgn_y * np.arccos(alpha*a[0] / np.linalg.norm(alpha * a[[0, 1]]))
    
    # 角マーク用の角bの値(球面座標)を計算
    psi_b = np.arccos(beta*b[2] / np.linalg.norm(beta * b))
    sgn_y = 1.0 if beta*b[1] >= 0.0 else -1.0
    phi_b = sgn_y * np.arccos(beta*b[0] / np.linalg.norm(beta * b[[0, 1]]))
    
    # 差が180°を超える場合は、負の角度を優角に変換
    if abs(phi_a - phi_b) > np.pi:
        phi_a = phi_a if alpha*a[1] >= 0.0 else phi_a + 2.0*np.pi
        phi_b = phi_b if beta*b[1] >= 0.0 else phi_b + 2.0*np.pi
    
    # 角マーク用の値(ラジアン)を作成
    rad_psi_vals = np.linspace(start=psi_a, stop=psi_b, num=100)
    rad_phi_vals = np.linspace(start=phi_a, stop=phi_b, num=100)

    # 角マーク用の座標を計算
    r = 0.2
    angle_x_vals = np.sin(rad_psi_vals) * np.cos(rad_phi_vals) * r
    angle_y_vals = np.sin(rad_psi_vals) * np.sin(rad_phi_vals) * r
    angle_z_vals = np.cos(rad_psi_vals) * r
    angle_mark_theta = (angle_x_vals, angle_y_vals, angle_z_vals)

    # 角ラベル用の座標を計算
    r = 0.4
    angle_x = np.sin(np.median(rad_psi_vals)) * np.cos(np.median(rad_phi_vals)) * r
    angle_y = np.sin(np.median(rad_psi_vals)) * np.sin(np.median(rad_phi_vals)) * r
    angle_z = np.cos(np.median(rad_psi_vals)) * r
    angle_label_theta = (angle_x, angle_y, angle_z)
    
    return angle_mark_theta, angle_label_theta, angle_mark_xy, angle_mark_xz

　詳しくは「2つのベクトルのなす角の可視化」を参照してください。

　3次元ベクトルを指定します。

# ベクトルを指定
a = np.array([2.0, 1.5, -1.0])
b = np.array([-1.5, 1.5, 1.0])

# 標準偏差を計算
std_a = np.std(a)
std_b = np.std(b)

# 平均除去ベクトルに変換
a -= np.mean(a)
b -= np.mean(b)
print(a)
print(b)

[ 1.16666667  0.66666667 -1.83333333]
[-1.83333333  1.16666667  0.66666667]

　 $\tilde{\mathbf{a}} = (\tilde{a}_1, \tilde{a}_2, \tilde{a}_3)^{\top}$ をa、 $\tilde{\mathbf{b}} = (\tilde{b}_1, \tilde{b}_2, \tilde{b}_3)^{\top}$ をbとして値を指定します。ただし、Pythonでは0からインデックスが割り当てられるので、 $\tilde{a}_1$ の値はa[0]に対応します。

　(平均が0のベクトルを指定するのは面倒なので)2つのベクトルの値を指定して、それぞれ平均 $\mathrm{avg}(\mathbf{x})$ を引きます。本の表記に合わせると、元のベクトルが $\mathbf{a}, \mathbf{b}$ で、平均除去ベクトルが $\tilde{\mathbf{a}}, \tilde{\mathbf{b}}$ で、次の関係で表されます。

$\displaystyle \begin{aligned} \tilde{\mathbf{a}} &= \mathbf{a} - \mathrm{avg}(\mathbf{a}) \mathbf{1} \\ \tilde{\mathbf{b}} &= \mathbf{b} - \mathrm{avg}(\mathbf{b}) \mathbf{1} \end{aligned}$

　元のベクトルは使わないので、プログラム上の変数名はどちらもa, bとします。
　また、標準化の処理に使うため、元のベクトルの標準偏差 $\mathrm{std}(\mathbf{x})$ を保存しておきます。

　2つのベクトルのなす角を計算します。

# 平均除去ベクトルa,bのなす角を計算
dot_ab = np.dot(a, b)
norm_a = np.linalg.norm(a)
norm_b = np.linalg.norm(b)
cos_theta = dot_ab / norm_a / norm_b
theta = np.arccos(cos_theta)
print(cos_theta) # 相関係数
print(theta) # ラジアン

-0.5
2.0943951023931957

　 $\tilde{\mathbf{a}}$ と $\tilde{\mathbf{b}}$ のなす角 $\theta$ は、次の式で定義されます。

$\displaystyle \begin{aligned} \cos \theta &= \frac{ \tilde{\mathbf{a}}^{\top} \tilde{\mathbf{b}} }{ \|\tilde{\mathbf{a}}\| \|\tilde{\mathbf{b}}\| } \\ \theta &= \arccos \left( \frac{ \tilde{\mathbf{a}}^{\top} \tilde{\mathbf{b}} }{ \|\tilde{\mathbf{a}}\| \|\tilde{\mathbf{b}}\| } \right) \end{aligned}$

　コサイン関数 $\cos x$ の逆関数(逆コサイン関数) $\arccos x$ はnp.arccos()、ユークリッドノルム $\|\mathbf{x}\|$ はnp.linalg.norm()で計算できます。
　ただし、 $\tilde{\mathbf{a}}$ または $\tilde{\mathbf{b}}$ が0ベクトルだと0除算になるため定義できません(計算結果がnanになります)。

　ベクトルを標準化して、相関係数を計算します。

# 標準化
u = a / std_a
v = b / std_b

# 相関係数を計算
rho = np.dot(u, v) / len(u)
print(rho) # 相関係数
print(np.arccos(rho)) # ラジアン

-0.5000000000000001
2.0943951023931957

　平均除去ベクトル(の各要素)を元のベクトルの標準偏差で割って、標準化します。

$\displaystyle \begin{aligned} \mathbf{u} &= \frac{ \tilde{\mathbf{a}} }{ \mathrm{std}(\mathbf{a}) } = \frac{ \mathbf{a} - \mathrm{avg}(\mathbf{a}) \mathbf{1} }{ \mathrm{std}(\mathbf{a}) } \\ \mathbf{v} &= \frac{ \tilde{\mathbf{b}} }{ \mathrm{std}(\mathbf{b}) } = \frac{ \mathbf{b} - \mathrm{avg}(\mathbf{b}) \mathbf{1} }{ \mathrm{std}(\mathbf{b}) } \end{aligned}$

　 $\mathbf{a}, \mathbf{b}$ を標準化したベクトル $\mathbf{u}, \mathbf{v}$ の内積を次元数 $D$ で割って、相関係数を計算します。

$\displaystyle \rho = \frac{\mathbf{u}^{\top} \mathbf{v}}{D}$

　標準化ベクトルを用いて計算した相関係数rhoと、平均除去ベクトルを用いて計算したcos_thetaが一致するのを確認できます。

　3次元空間上の「2つのベクトル $\tilde{\mathbf{a}}, \tilde{\mathbf{b}}$ のなす角 $\theta$ 」のグラフ、コサイン曲線上の「なす角 $\theta$ と相関係数 $\rho = \cos \theta$ 」のグラフ、「ベクトル $\tilde{\mathbf{a}}, \tilde{\mathbf{b}}$ の成分」のグラフ、「 $\tilde{\mathbf{a}}, \tilde{\mathbf{b}}$ の散布図と相関係数の直線 $y = \rho x$ 」のグラフを並べて作成します。

・作図コード(クリックで展開)

# 3Dベクトルのなす角を作図
fig = plt.figure(figsize=(12, 10), facecolor='white')
fig.suptitle('angle and correlation coefficient', fontsize=20)

# なす角の作図用の配列を作成
angle_mark_theta, angle_label_theta, angle_mark_xy, angle_mark_xz = angle3d(a, b)
angle_xy_x_vals, angle_xy_y_vals = angle_mark_xy
angle_xz_x_vals, angle_xz_z_vals = angle_mark_xz
angle_x_vals, angle_y_vals, angle_z_vals = angle_mark_theta
angle_x, angle_y, angle_z = angle_label_theta

# グラフサイズ用の値を設定
x_min = np.floor(np.min([0.0, a[0], b[0]])) - 0.1
x_max = np.ceil(np.max([0.0, a[0], b[0]])) + 0.1
y_min = np.floor(np.min([0.0, a[1], b[1]])) - 0.1
y_max = np.ceil(np.max([0.0, a[1], b[1]])) + 0.1
z_min = np.floor(np.min([0.0, a[2], b[2]])) - 0.1
z_max = np.ceil(np.max([0.0, a[2], b[2]])) + 0.1

# 3Dベクトルのなす角を作図
ax00 = fig.add_subplot(2, 2, 1, projection='3d') # 左上
ax00.quiver(0, 0, 0, *a, color='red', arrow_length_ratio=0.1, 
            label='$\\tilde{a} = ('+', '.join(map(str, a.round(2)))+')$') # ベクトルa
ax00.quiver(0, 0, 0, *b, color='blue', arrow_length_ratio=0.1, 
            label='$\\tilde{b} = ('+', '.join(map(str, b.round(2)))+')$') # ベクトルb
ax00.quiver(0, y_max, 0, a[0], 0, a[2], color='red', linestyle='--', arrow_length_ratio=0.0) # xy面上のベクトルa
ax00.quiver(0, 0, z_min, a[0], a[1], 0, color='red', linestyle='--', arrow_length_ratio=0.0) # xz面上のベクトルa
ax00.quiver(0, y_max, 0, b[0], 0, b[2], color='blue', linestyle='--', arrow_length_ratio=0.0) # xy面上のベクトルb
ax00.quiver(0, 0, z_min, b[0], b[1], 0, color='blue', linestyle='--', arrow_length_ratio=0.0) # xz面上のベクトルb
ax00.plot(angle_xy_x_vals, angle_xy_y_vals, np.repeat(z_min, len(angle_xy_x_vals)), 
          color='black', linewidth=1) # xy面上の角マーク
ax00.plot(angle_xz_x_vals, np.repeat(y_max, len(angle_xz_x_vals)), angle_xz_z_vals, 
          color='black', linewidth=1) # xz面上の角マーク
ax00.plot(angle_x_vals, angle_y_vals, angle_z_vals, color='black', linewidth=1) # 角マーク
ax00.text(angle_x, angle_y, angle_z, s='$\\theta$', size=15, ha='center', va='center') # 角ラベル
ax00.quiver([0, a[0], b[0], 0, a[0], b[0]], 
            [0, a[1], b[1], y_max, y_max, y_max], 
            [z_min, z_min, z_min, 0, a[2], b[2]], 
            [0, 0, 0, 0, 0, 0], 
            [0, 0, 0, -y_max, a[1]-y_max, b[1]-y_max], 
            [-z_min, a[2]-z_min, b[2]-z_min, 0, 0, 0], 
            color='gray', arrow_length_ratio=0, linestyle=':') # 補助線
ax00.set_xlim(x_min, x_max)
ax00.set_ylim(y_min, y_max)
ax00.set_zlim(z_min, z_max)
ax00.set_xlabel('$x_1$')
ax00.set_ylabel('$x_2$')
ax00.set_zlabel('$x_3$')
ax00.set_title('$\\theta=' + str(theta.round(2)) + ', ' + 
               '\\theta^\circ=' + str((theta * 180.0 / np.pi).round(1)) + '$', loc='left')
ax00.legend()
ax00.set_aspect('equal')
#ax00.view_init(elev=90, azim=270) # xy面
#ax00.view_init(elev=0, azim=270) # xz面

# コサイン曲線を計算
theta_vals = np.linspace(start=0.0, stop=2.0*np.pi, num=100)
cos_vals = np.cos(theta_vals)

# 軸ラベル用の値を設定
denom = 3.0
n_ticks = np.arange(start=0, stop=2.0*denom+1, step=1)
t_labels = ['$\\frac{'+str(int(n))+'}{'+str(int(denom))+'} \pi$' for n in n_ticks]

# コサイン曲線(なす角と相関係数の関係)を作図
ax01 = fig.add_subplot(2, 2, 2) # 右上
ax01.plot(theta_vals, cos_vals) # コサイン曲線
ax01.scatter(theta, cos_theta, color='black', s=100, label='$(\\theta, \cos \\theta)$') # 曲線上の点
ax01.set_xticks(ticks=n_ticks*np.pi/denom, labels=t_labels)
ax01.grid()
ax01.set_xlabel('$\\theta$')
ax01.set_ylabel('$\cos \\theta$')
ax01.set_title('$\\theta=' + str((theta / np.pi).round(2)) + '\pi, ' + 
               '\cos \\theta=' + str(cos_theta.round(2)) + '$', loc='left')
ax01.legend()

# x軸の値を作成
d_vals = np.arange(len(a))

# グラフサイズ用の値を設定
axis_size = np.max([*abs(a), *abs(b)]) + 0.5

# 平均除去ベクトルの各成分を作図
ax10 = fig.add_subplot(2, 2, 3) # 左下
ax10.plot(d_vals, a, color='red') # ベクトルaの成分の線
ax10.plot(d_vals, b, color='blue') # ベクトルbの成分の線
ax10.scatter(d_vals, a, color='red', s=100, label='$(d, \\tilde{a}_d)$') # ベクトルaの成分の点
ax10.scatter(d_vals, b, color='blue', s=100, label='$(d, \\tilde{b}_d)$') # ベクトルbの成分の点
ax10.set_xticks(ticks=d_vals, labels=d_vals+1)
ax10.set_ylim(bottom=-axis_size, top=axis_size)
ax10.grid()
ax10.set_xlabel('d')
ax10.set_ylabel('$\\tilde{a}_d, \\tilde{b}_d$')
ax10.legend()

# 相関係数の直線を作成
x_vals = np.linspace(start=-axis_size, stop=axis_size, num=100)
y_vals = rho * x_vals

# 平均除去ベクトルの散布図と相関係数を作図
ax11 = fig.add_subplot(2, 2, 4) # 右下
ax11.scatter(a, b, color='purple', s=100, label='$(\\tilde{a}_d, \\tilde{b}_d)$') # 次元ごとの成分の点
ax11.plot(x_vals, y_vals, color='black', label='$y = \\rho x$') # 相関係数の直線
ax11.set_xlabel('$\\tilde{a}_d$')
ax11.set_ylabel('$\\tilde{b}_d$')
ax11.set_title('$\\rho=' + str(rho.round(2)) + '$', loc='left')
ax11.set_xlim(left=-axis_size, right=axis_size)
ax11.set_ylim(bottom=-axis_size, top=axis_size)
ax11.grid()
ax11.legend()
ax11.set_aspect('equal')

plt.show()

　左上の図(1番目のサブプロット)については「2つのベクトルのなす角の可視化」を参照してください。

　右上の図(2番目のサブプロット)は、x軸を $0 \leq t \leq 2 \pi$ のラジアン(度数法の角度で表すと $0^{\circ} \leq t^{\circ} \leq 360^{\circ}$ )、y軸を $\cos t$ とするコサイン曲線を描画します。
　また曲線上に、なす角 $\theta$ に関する点を描画します。y軸の値は、cos_thetaまたはrhoを使います。

　左下の図(3番目のサブプロット)は、x軸を次元番号 $d$ 、y軸を $\tilde{a}_d, \tilde{b}_d$ として、散布図と折れ線グラフを描画します。

　右下の図(4番目のサブプロット)は、x軸を $\tilde{a}_d$ 、y軸を $\tilde{b}_d$ として散布図を描画します。
　また、相関係数 $\rho$ を傾きとした直線 $y = \rho x$ を描画します。

　ベクトルの値を変化させたアニメーションを作成します。

・作図コード(クリックで展開)

# フレーム数を設定
frame_num = 150

# ベクトルとして利用する値を指定
psi_a_vals = np.linspace(start=-4.0/3.0*np.pi, stop=8.0/3.0*np.pi, num=frame_num+1)[:frame_num]
phi_a_vals = np.linspace(start=5.0/2.0*np.pi, stop=-7.0/2.0*np.pi, num=frame_num+1)[:frame_num]
a_vals = np.array(
    [np.sin(psi_a_vals) * np.cos(phi_a_vals), 
     np.sin(psi_a_vals) * np.sin(phi_a_vals), 
     np.cos(psi_a_vals)]
).T
psi_b_vals = np.linspace(start=-9.0/4.0*np.pi, stop=7.0/4.0*np.pi, num=frame_num+1)[:frame_num]
phi_b_vals = np.linspace(start=0.0/5.0*np.pi, stop=10.0/5.0*np.pi, num=frame_num+1)[:frame_num]
b_vals = np.array(
    [np.sin(psi_b_vals) * np.cos(phi_b_vals), 
     np.sin(psi_b_vals) * np.sin(phi_b_vals), 
     np.cos(psi_b_vals)]
).T

# 平均除去ベクトルに変換
a_vals -= np.mean(a_vals, axis=1, keepdims=True)
b_vals -= np.mean(b_vals, axis=1, keepdims=True)

# グラフサイズ用の値を設定
x_min = np.floor(np.min([0.0, *a_vals[:, 0], *b_vals[:, 0]])) - 0.1
x_max = np.ceil(np.max([0.0, *a_vals[:, 0], *b_vals[:, 0]])) + 0.1
y_min = np.floor(np.min([0.0, *a_vals[:, 1], *b_vals[:, 1]])) - 0.1
y_max = np.ceil(np.max([0.0, *a_vals[:, 1], *b_vals[:, 1]])) + 0.1
z_min = np.floor(np.min([0.0, *a_vals[:, 2], *b_vals[:, 2]])) - 0.1
z_max = np.ceil(np.max([0.0, *a_vals[:, 2], *b_vals[:, 2]])) + 0.1
u_max = np.ceil((abs(a_vals) / np.std(a_vals, axis=1, keepdims=True)).max())
v_max = np.ceil((abs(b_vals) / np.std(b_vals, axis=1, keepdims=True)).max())
axis_size = np.max([u_max, v_max]) + 0.1

# コサイン波を計算
theta_vals = np.linspace(start=0.0, stop=2.0*np.pi, num=100)
cos_vals = np.cos(theta_vals)

# 軸ラベル用の値を設定
denom = 3.0
n_ticks = np.arange(start=0, stop=2.0*denom+1, step=1)
t_labels = ['$\\frac{'+str(int(n))+'}{'+str(int(denom))+'} \pi$' for n in n_ticks]
d_vals = np.arange(len(a_vals[0]))

# 作図用のオブジェクトを初期化
fig = plt.figure(figsize=(12, 10), facecolor='white')
fig.suptitle('angle and correlation coefficient', fontsize=20)
ax00 = fig.add_subplot(2, 2, 1, projection='3d') # 左上
ax01 = fig.add_subplot(2, 2, 2) # 右上
ax10 = fig.add_subplot(2, 2, 3) # 左下
ax11 = fig.add_subplot(2, 2, 4) # 右下

# 作図処理を関数として定義
def update(i):
    
    # 前フレームのグラフを初期化
    ax00.clear()
    ax01.clear()
    ax10.clear()
    ax11.clear()
    
    # i番目の平均除去ベクトルを取得
    a = a_vals[i]
    b = b_vals[i]
    
    # 相関係数を計算
    rho = np.dot(a, b) / np.linalg.norm(a) / np.linalg.norm(b)
    
    # 平均除去ベクトルa,bのなす角を計算
    theta = np.arccos(rho)
    
    # なす角の作図用の配列を作成
    angle_mark_theta, angle_label_theta, angle_mark_xy, angle_mark_xz = angle3d(a, b)
    angle_xy_x_vals, angle_xy_y_vals = angle_mark_xy
    angle_xz_x_vals, angle_xz_z_vals = angle_mark_xz
    angle_x_vals, angle_y_vals, angle_z_vals = angle_mark_theta
    angle_x, angle_y, angle_z = angle_label_theta
    
    # 3Dベクトルのなす角を作図
    ax00.quiver(0, 0, 0, *a, color='red', arrow_length_ratio=0.1, 
                label='$\\tilde{a} = ('+', '.join(map(str, a.round(2)))+')$') # ベクトルa
    ax00.quiver(0, 0, 0, *b, color='blue', arrow_length_ratio=0.1, 
                label='$\\tilde{b} = ('+', '.join(map(str, b.round(2)))+')$') # ベクトルb
    ax00.quiver(0, y_max, 0, a[0], 0, a[2], color='red', linestyle='--', arrow_length_ratio=0.0) # xy面上のベクトルa
    ax00.quiver(0, 0, z_min, a[0], a[1], 0, color='red', linestyle='--', arrow_length_ratio=0.0) # xz面上のベクトルa
    ax00.quiver(0, y_max, 0, b[0], 0, b[2], color='blue', linestyle='--', arrow_length_ratio=0.0) # xy面上のベクトルb
    ax00.quiver(0, 0, z_min, b[0], b[1], 0, color='blue', linestyle='--', arrow_length_ratio=0.0) # xz面上のベクトルb
    ax00.plot(angle_xy_x_vals, angle_xy_y_vals, np.repeat(z_min, len(angle_xy_x_vals)), 
              color='black', linewidth=1) # xy面上の角マーク
    ax00.plot(angle_xz_x_vals, np.repeat(y_max, len(angle_xz_x_vals)), angle_xz_z_vals, 
              color='black', linewidth=1) # xz面上の角マーク
    ax00.plot(angle_x_vals, angle_y_vals, angle_z_vals, color='black', linewidth=1) # 角マーク
    ax00.text(angle_x, angle_y, angle_z, s='$\\theta$', size=15, ha='center', va='center') # 角ラベル
    ax00.quiver([0, a[0], b[0], 0, a[0], b[0]], 
                [0, a[1], b[1], y_max, y_max, y_max], 
                [z_min, z_min, z_min, 0, a[2], b[2]], 
                [0, 0, 0, 0, 0, 0], 
                [0, 0, 0, -y_max, a[1]-y_max, b[1]-y_max], 
                [-z_min, a[2]-z_min, b[2]-z_min, 0, 0, 0], 
                color='gray', arrow_length_ratio=0, linestyle=':') # 補助線
    ax00.set_xlim(x_min, x_max)
    ax00.set_ylim(y_min, y_max)
    ax00.set_zlim(z_min, z_max)
    ax00.set_xlabel('$x_1$')
    ax00.set_ylabel('$x_2$')
    ax00.set_zlabel('$x_3$')
    ax00.set_title('$\\theta=' + str(theta.round(2)) + ', ' + 
                   '\\theta^\circ=' + str((theta * 180.0 / np.pi).round(1)) + '$', loc='left')
    ax00.legend(loc='upper left')
    ax00.set_aspect('equal')
    #ax00.view_init(elev=90, azim=270) # xy面
    #ax00.view_init(elev=0, azim=270) # xz面

    # コサイン曲線(なす角と相関係数の関係)を作図
    ax01.plot(theta_vals, cos_vals) # コサイン曲線
    ax01.scatter(theta, rho, color='black', s=100, label='$(\\theta, \cos \\theta)$') # 曲線上の点
    ax01.set_xticks(ticks=n_ticks*np.pi/denom, labels=t_labels)
    ax01.grid()
    ax01.set_xlabel('$\\theta$')
    ax01.set_ylabel('$\cos \\theta$')
    ax01.set_title('$\\theta=' + str((theta / np.pi).round(2)) + '\pi, ' + 
                   '\cos \\theta=' + str(rho.round(2)) + '$', loc='left')
    ax01.legend(loc='lower left')
    
    # 平均除去ベクトルの各成分を作図
    ax10.plot(d_vals, a, color='red') # ベクトルaの成分の線
    ax10.plot(d_vals, b, color='blue') # ベクトルbの成分の線
    ax10.scatter(d_vals, a, color='red', s=100, label='$(d, \\tilde{a}_d)$') # ベクトルaの成分の点
    ax10.scatter(d_vals, b, color='blue', s=100, label='$(d, \\tilde{b}_d)$') # ベクトルbの成分の点
    ax10.set_xticks(ticks=d_vals, labels=d_vals+1)
    ax10.set_ylim(bottom=-axis_size, top=axis_size)
    ax10.grid()
    ax10.set_xlabel('d')
    ax10.set_ylabel('$\\tilde{a}_d, \\tilde{b}_d$')
    ax10.legend(loc='upper left')

    # 相関係数の直線を作成
    x_vals = np.linspace(start=-axis_size, stop=axis_size, num=100)
    y_vals = rho * x_vals

    # 平均除去ベクトルの散布図と相関係数を作図
    ax11.scatter(a, b, color='purple', s=100, label='$(\\tilde{a}_d, \\tilde{b}_d)$') # 次元ごとの成分の点
    ax11.plot(x_vals, y_vals, color='black', label='$y = \\rho x$') # 相関係数の直線
    ax11.set_xlim(left=-axis_size, right=axis_size)
    ax11.set_ylim(bottom=-axis_size, top=axis_size)
    ax11.grid()
    ax11.set_xlabel('$\\tilde{a}_d$')
    ax11.set_ylabel('$\\tilde{b}_d$')
    ax11.set_title('$\\rho=' + str(rho.round(2)) + '$', loc='left')
    ax11.legend(loc='upper left')    
    ax11.set_aspect('equal')

# gif画像を作成
ani = FuncAnimation(fig=fig, func=update, frames=frame_num, interval=100)

# gif画像を保存
ani.save('corr_3d.gif')

　作図処理をupdate()として定義して、FuncAnimation()でgif画像を作成します。

　なす角は $0 \leq \theta \leq \pi$ の範囲(度数法で表すと $0^{\circ} \leq \theta^{\circ} \leq 180^{\circ}$ )の値をとるので、コサイン曲線の左半分を点が変化します。ちなみに、なす角と反対側の角度に注目すると、 $2 \pi - \theta$ (度数法だと $360^{\circ} - \theta^{\circ}$ )であり、 $\pi$ から $2 \pi$ の値をとります(右半分の点になります)。コサイン曲線は $x = \pi$ で対称な形なので、 $\cos \theta = \cos (2 \pi - \theta)$ です。つまり、どちらの角度に注目しても相関係数が同じ値になるのが分かります。
　コサイン関数は $-1 \leq \cos \theta \leq 1$ の値をとるので、相関係数も $-1 \leq \rho \leq 1$ の範囲になります。
　また、2つのベクトルの向き(折れ線の形)が近いほど、なす角が小さく、相関係数が大きくなり(1に近付き)ます。逆に2つのベクトルの向きが遠い(折れ線の形が異なる)ほど、なす角が大きく、相関係数が小さくなり(-1に近付き)ます。

　ちなみに、点 $\tilde{\mathbf{a}}, \tilde{\mathbf{b}}$ (矢印の先端)は、次の軌道をしています。

# ベクトルの軌道を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.plot(*a_vals.T, color='red', alpha=0.5, label='$\\tilde{a}$') # ベクトルaの軌道
ax.plot(*b_vals.T, color='blue', alpha=0.5, label='$\\tilde{b}$') # ベクトルbの軌道
ax.set_xticks(ticks=np.arange(x_min, x_max+1))
ax.set_yticks(ticks=np.arange(y_min, y_max+1))
ax.set_zticks(ticks=np.arange(z_min, z_max+1))
ax.set_ylim(y_min, y_max)
ax.set_zlim(z_min, z_max)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
#ax.set_title('$\\beta='+str(beta)+'$', loc='left')
fig.suptitle('trace $\\tilde{a}, \\tilde{b}$', fontsize=20)
ax.legend()
ax.set_aspect('equal')
plt.show()

　a_vals, b_valsを転置して、次元ごとにアンパックして描画します。
(対角線上に伸びる形？になる理由が分かりません。)

　この記事では、なす角と相関係数を可視化しました。次の記事では、k平均法によるクラスタリングを実装します。

参考書籍

Stephen Boyd・Lieven Vandenberghe(著),玉木徹(訳)『スタンフォード　ベクトル・行列からはじめる最適化数学』講談社サイエンティク,2021年.

おわりに

　1章で扱ったものの、どうしてもベクトルというものを感覚的に掴めていなかった気がします。相関係数の説明でよく見る散布図に関しても、対角線上に散らばるn個の点をn個のデータと思っていました。各点が1つのデータというのもそうなんでしょうが、x軸とy軸の値がそれぞれ別のn次元ベクトルだったとは、そんな風には考えたことはありませんでした。
　そんなこんなを自分の頭に刷り込みたいと思ってこの記事を書き(図を作り)ました。

　上記のベクトル表現の理解のために、3次元空間上に(文字通りの)ベクトルを2つ描画して、その座標を平面上に3点描画しました。3点だと散布図っぽくなかったですがしょうがない。
　それとは別にベクトルの類似度を、ベクトルの向き・なす角の大小・折れ線グラフの形状・散布図の形で表現しました。折れ線が山型と谷型のように大小が逆の形だと左上から右下に点が並び、単にバラバラだと散布図もバラバラになり、同じ形だと左下から右上に並びます。改めて考えるとそりゃそうなんですが、新鮮でした。
　さらに相関係数が-1から1になるのが分かるように、なす角が $\pi$ 以上にならないのを示しつつ、コサイン曲線上の点を描画しました。
　以上の4つの図で、頭の中でふわっとしていたアレコレを理解できました。

　しかし、理解の過程で(図が完成するまでに)色々試行錯誤しました。
　まず、ベクトルから平均を引くとどうなるのかも分かっていませんでした。向きも変わるんですね(正負も変わったりするんだから当然)。というわけで、平均除去ベクトルから話がスタートします。そして、散布図には平均除去ベクトルなのか標準化ベクトルなのか、どのベクトルの値を使えばいいのかも悩みました。傾きが相関係数で原点を通る直線と対応させたかったので、中心が原点になるような値になればいいので、平均除去ベクトルでも標準化ベクトルでも問題なさそうです。他の図との対応から平均除去ベクトルを使いました。
　他にも、1周期分のコサイン曲線を描画するためにx軸の範囲を $2 \pi$ にしたら、点が左半分しか移動せずなんでだとなり、いやなす角は $\pi$ までだわと気付き、このあとがきを書いてたら反対側の角度が右半分なのかと気付きました。いやぁ手癖で1周期分描いといてよかった。
　手癖になった理由はこのシリーズ「Rによる三角関数(円関数)入門：記事一覧 - からっぽのしょこ」です。こっちも、もっとちゃんと勉強しとけばよかったシリーズなので、合わせてどうですか。

　とまぁアレコレ頭を捻ったこともあって、このシリーズのこれまでの記事で一番刺激的で面白く満足できた内容になりました。結構苦労した3次元なす角の可視化も再利用できて良かったです。
　他の人に伝わるのかはまた別ですがね(笑)

　2023年3月20日は、エビ中こと私立恵比寿中学の小久保柚乃さんの16歳のお誕生日です！

　ここ半年ほどで完全に嵌りました。とても良いグループ・楽曲なので、ハロプロと共にぜひ聴いてみてください。

　きっかけは、「関内デビル」にハロメンが出るということでゲスト回の前後も含めて観てたら、その中(たしか司会の方がしゃべってるのに小道具のおもちゃを触ってたら電気ショックのやつだったゆのぴ回)でどタイプのメンバーがいるなぁとエビ中を調べてみて、「Anytime, Anywhere」を聴いてこれ私ずっと前から好きでしたと言うくらい良くって、なんで今まで知らなかったんだよ自分と思ったら新メン・新曲だと知り、少しずつ追うようになって、ハロプログループも参加した合同ライブで生パフォーマンスを観て完全に嵌りました。
　特にゆのぴのマイペースで不思議な雰囲気を醸し出しているけどいわゆる天然さんや不思議ちゃんではない感じとあと可愛いお顔と、エビ中の(ハロプロとはまた違った)楽し気な雰囲気を作りつつ地力がちゃんと凄いところが好きです。(まだそれほど知らないんですが。)

　というわけで今日からこのブログは、ハロオタときどきエビオタによる学習ブログとなります。今後ともよろしくお願いします。

【次の内容】