はじめに
素直なやり方ではできなかったのでむりくりなんとかする黒魔術シリーズです。もっといい方法があれば教えてください。
この記事では、Pythonで3次元の三角図を作成します。
【前の内容】
【目次】
Matplotlibで3D三角グラフを作図したい
前回は、三角座標上に等高線図(2次元のグラフ)を描画した。今回は、三角図上に曲面図(3次元のグラフ)を描画する。
三角図の座標の作図については前々回の記事、作図用と計算用の値の作成については前回の記事を参照のこと。
利用するライブラリを読み込む。
# 利用ライブラリ import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation from scipy.stats import dirichlet
3D三角座標
まずは、3次元三角図の基となる軸やグリッド線を用意する。
三角座標を描画する用の配列を作成する。詳しくは「Matplotlibで三角グラフを作図したい - からっぽのしょこ」を参照のこと。
# 軸目盛の位置を指定 axis_vals = np.arange(start=0.0, stop=1.1, step=0.1) # 軸線用の値を作成 axis_x = np.array([0.5, 0.0, 1.0]) axis_y = np.array([0.5*np.sqrt(3.0), 0.0, 0.0]) axis_u = np.array([-0.5, 1.0, -0.5]) axis_v = np.array([-0.5*np.sqrt(3.0), 0.0, 0.5*np.sqrt(3.0)]) # グリッド線用の値を作成 grid_x = np.hstack([ 0.5 * axis_vals, axis_vals, 0.5 * axis_vals + 0.5 ]) grid_y = np.hstack([ 0.5 * axis_vals * np.sqrt(3.0), np.zeros_like(axis_vals), 0.5 * (1.0 - axis_vals) * np.sqrt(3.0) ]) grid_u = np.hstack([ 0.5 * axis_vals, 0.5 * (1.0 - axis_vals), -axis_vals ]) grid_v = np.hstack([ -0.5 * axis_vals * np.sqrt(3.0), 0.5 * (1.0 - axis_vals) * np.sqrt(3.0), np.zeros_like(axis_vals) ])
作図用の格子点と計算用の3次元変数の値を作成する。詳しくは「Matplotlibで三角グラフの等高線を作図したい - からっぽのしょこ」を参照のこと。
# 2次元座標の値を作成 y_0_vals = np.linspace(start=0.0, stop=1.0, num=201) y_1_vals = np.linspace(start=0.0, stop=0.5*np.sqrt(3.0), num=200) # 2次元座標の格子点を作成 y_0_grid, y_1_grid = np.meshgrid(y_0_vals, y_1_vals) # 格子点の形状を保存 y_shape = y_0_grid.shape # 3次元座標の値に変換 x_1_vals = y_0_grid.flatten() - y_1_grid.flatten() / np.sqrt(3.0) x_2_vals = 2.0 * y_1_grid.flatten() / np.sqrt(3.0) # 範囲外の点を欠損値に置換 x_1_vals = np.where( (x_1_vals >= 0.0) & (x_1_vals <= 1.0), x_1_vals, np.nan ) x_2_vals = np.where( (x_2_vals >= 0.0) & (x_2_vals <= 1.0), x_2_vals, np.nan ) # 3次元座標の値に変換 x_0_vals = 1.0 - x_1_vals - x_2_vals # 範囲外の点を欠損値に置換 x_0_vals = np.where( (x_0_vals >= 0.0) & (x_0_vals <= 1.0), x_0_vals, np.nan ) # 計算用の3次元座標の点を作成 x_points = np.stack([x_0_vals, x_1_vals, x_2_vals], axis=1)
3次元の三角図の基となる軸やグリッド線を用意する。
# 三角座標上の曲面図を作成 fig = plt.figure(figsize=(12, 10), facecolor='white') # 図の設定 ax = fig.add_subplot(projection='3d') # 3D用の設定 ax.quiver(grid_x, grid_y, np.zeros_like(grid_x), grid_u, grid_v, np.zeros_like(grid_x), arrow_length_ratio=0.0, ec='gray', linewidth=1.5, linestyle=':') # 三角座標のグリッド線 ax.quiver(axis_x, axis_y, np.zeros_like(axis_x), axis_u, axis_v, np.zeros_like(axis_x), arrow_length_ratio=0.0, ec='black', linestyle='-') # 三角座標の枠線 for val in axis_vals: ax.text(x=0.5*val-0.05, y=0.5*val*np.sqrt(3.0), z=0.0, s=str(np.round(1.0-val, 1)), ha='center', va='center') # 三角座標のx軸目盛 ax.text(x=val, y=0.0-0.05, z=0.0, s=str(np.round(val, 1)), ha='center', va='center') # 三角座標のy軸目盛 ax.text(x=0.5*val+0.5+0.05, y=0.5*(1.0-val)*np.sqrt(3.0), z=0.0, s=str(np.round(1.0-val, 1)), ha='center', va='center') # 三角座標のz軸目盛 ax.text(x=0.25-0.1, y=0.25*np.sqrt(3.0), z=0.0, s='$x_0$', ha='right', va='center', size=25) # 三角座標のx軸ラベル ax.text(x=0.5, y=0.0-0.1, z=0.0-0.1, s='$x_1$', ha='center', va='top', size=25) # 三角座標のy軸ラベル ax.text(x=0.75+0.1, y=0.25*np.sqrt(3.0), z=0.0, s='$x_2$', ha='left', va='center', size=25) # 三角図のz軸ラベル ax.set_xticks(ticks=[0.0, 0.5, 1.0], labels='') # 2次元座標のx軸目盛 ax.set_yticks(ticks=[0.0, 0.25*np.sqrt(3.0), 0.5*np.sqrt(3.0)], labels='') # 2次元座標のy軸目盛 ax.set_zlabel(zlabel='$f(x_0, x_1, x_2)$') # z軸ラベル ax.set_zlim(bottom=0.0, top=1.0) # z軸の表示範囲 ax.set_box_aspect(aspect=(1, 1, 1)) # アスペクト比 fig.suptitle(t='3D Ternary Plot', fontsize=20) # タイトル #ax.view_init(elev=90, azim=-90) # 表示角度 plt.show() # 描画
次からは、ここに曲面図を描画する。
曲面図
三角座標上の曲面図を作成する。例として、ディリクレ分布の確率密度の曲面図を作成する。
ディリクレ分布のパラメータを設定して、確率密度を計算する。
# ディリクレ分布のパラメータを指定 alpha_k = np.array([2.0, 4.0, 6.0]) # ディリクレ分布の確率密度を計算 dens_vals = np.array( [dirichlet.pdf(x=x_k, alpha=alpha_k) if all(x_k != np.nan) else np.nan for x_k in x_points] ) print(dens_vals)
[ 0. 0. 0. ... nan nan nan]
ディリクレ分布の確率密度は、SciPyライブラリのディリクレ分布のモジュールdirichletのpdf()メソッドで計算できる。確率変数の引数xにx_pointsの各行、パラメータの引数alphaに設定したパラメータalpha_kを指定する。
リスト内包表記を使って、x_pointsの行ごとに確率密度を計算する。ただし、np.nanを含む行の場合は計算を行わず、np.nanを格納する。
三角図上に確率密度の曲面を描画する。
# パラメータラベル用の文字列を作成 param_text = '$' + '\\alpha=('+', '.join([str(val) for val in alpha_k])+')' + ', x=(x_0, x_1, x_2)' + '$' # 三角座標上の曲面図を作成 fig = plt.figure(figsize=(12, 10), facecolor='white') # 図の設定 ax = fig.add_subplot(projection='3d') # 3D用の設定 ax.quiver(grid_x, grid_y, np.zeros_like(grid_x), grid_u, grid_v, np.zeros_like(grid_x), arrow_length_ratio=0.0, ec='gray', linewidth=1.5, linestyle=':') # 三角座標のグリッド線 ax.quiver(axis_x, axis_y, np.zeros_like(axis_x), axis_u, axis_v, np.zeros_like(axis_x), arrow_length_ratio=0.0, ec='black', linestyle='-') # 三角座標の枠線 for val in axis_vals: ax.text(x=0.5*val-0.05, y=0.5*val*np.sqrt(3.0), z=0.0, s=str(np.round(1.0-val, 1)), ha='center', va='center') # 三角座標のx軸目盛 ax.text(x=val, y=0.0-0.05, z=0.0, s=str(np.round(val, 1)), ha='center', va='center') # 三角座標のy軸目盛 ax.text(x=0.5*val+0.5+0.05, y=0.5*(1.0-val)*np.sqrt(3.0), z=0.0, s=str(np.round(1.0-val, 1)), ha='center', va='center') # 三角座標のz軸目盛 ax.text(x=0.25-0.1, y=0.25*np.sqrt(3.0), z=0.0, s='$x_0$', ha='right', va='center', size=25) # 三角座標のx軸ラベル ax.text(x=0.5, y=0.0-0.1, z=0.0-0.1, s='$x_1$', ha='center', va='top', size=25) # 三角座標のy軸ラベル ax.text(x=0.75+0.1, y=0.25*np.sqrt(3.0), z=0.0, s='$x_2$', ha='left', va='center', size=25) # 三角図のz軸ラベル ax.contour(y_0_grid, y_1_grid, dens_vals.reshape(y_shape), offset=0.0) # 確率密度の等高線 ax.plot_surface(y_0_grid, y_1_grid, dens_vals.reshape(y_shape), cmap='viridis', alpha=0.8) # 確率密度の曲面 ax.set_xticks(ticks=[0.0, 0.5, 1.0], labels='') # 2次元座標のx軸目盛 ax.set_yticks(ticks=[0.0, 0.25*np.sqrt(3.0), 0.5*np.sqrt(3.0)], labels='') # 2次元座標のy軸目盛 ax.set_zlabel(zlabel='density') # z軸ラベル ax.set_box_aspect(aspect=(1, 1, 1)) # アスペクト比 fig.suptitle(t='3D Ternary Plot', fontsize=20) # タイトル ax.set_title(label=param_text, loc='left') # パラメータラベル #ax.view_init(elev=90, azim=-90) # 表示角度 plt.show() # 描画
ax.plot_surface()で曲面図を描画する。三角図外の点は欠損値np.nanなので描画されない(つまりデータの半分を捨てる非効率設計である)。
また、座標上に曲面の等高線も描画する。ax.contour()のoffset引数に表示するz軸の値を指定する。
以上で、私が欲しいグラフが得られた。
曲面図のアニメーション
最後に、等高線図のアニメーション(gif画像)を作成する。
水平方向に回転させてグラフを確認する。
# 水平方向の角度として利用する値を指定 h_vals = np.arange(0.0, 360.0, step=5.0) # フレーム数を設定 frame_num = len(h_vals) # パラメータラベル用の文字列を作成 param_text = '$' + '\\alpha=('+', '.join([str(val) for val in alpha_k])+')' + ', x=(x_0, x_1, x_2)' + '$' # 三角座標上の等高線図を作成 fig = plt.figure(figsize=(12, 12), facecolor='white') # 図の設定 ax = fig.add_subplot(projection='3d') # 3D用の設定 fig.suptitle(t='3D Ternary Plot', fontsize=20) # タイトル # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目の角度を取得 h = h_vals[i] # 三角座標上の曲面図を作成 ax.quiver(grid_x, grid_y, np.zeros_like(grid_x), grid_u, grid_v, np.zeros_like(grid_x), arrow_length_ratio=0.0, ec='gray', linewidth=1.5, linestyle=':') # 三角座標のグリッド線 ax.quiver(axis_x, axis_y, np.zeros_like(axis_x), axis_u, axis_v, np.zeros_like(axis_x), arrow_length_ratio=0.0, ec='black', linestyle='-') # 三角座標の枠線 for val in axis_vals: ax.text(x=0.5*val-0.05, y=0.5*val*np.sqrt(3.0), z=0.0, s=str(np.round(1.0-val, 1)), ha='center', va='center') # 三角座標のx軸目盛 ax.text(x=val, y=0.0-0.05, z=0.0, s=str(np.round(val, 1)), ha='center', va='center') # 三角座標のy軸目盛 ax.text(x=0.5*val+0.5+0.05, y=0.5*(1.0-val)*np.sqrt(3.0), z=0.0, s=str(np.round(1.0-val, 1)), ha='center', va='center') # 三角座標のz軸目盛 ax.text(x=0.25-0.1, y=0.25*np.sqrt(3.0), z=0.0, s='$x_0$', ha='right', va='center', size=25) # 三角座標のx軸ラベル ax.text(x=0.5, y=0.0-0.1, z=0.0-0.1, s='$x_1$', ha='center', va='top', size=25) # 三角座標のy軸ラベル ax.text(x=0.75+0.1, y=0.25*np.sqrt(3.0), z=0.0, s='$x_2$', ha='left', va='center', size=25) # 三角図のz軸ラベル ax.contour(y_0_grid, y_1_grid, dens_vals.reshape(y_shape), offset=0.0) # 確率密度の等高線 ax.plot_surface(y_0_grid, y_1_grid, dens_vals.reshape(y_shape), cmap='viridis', alpha=0.8) # 確率密度の曲面 ax.set_xticks(ticks=[0.0, 0.5, 1.0], labels='') # 2次元座標のx軸目盛 ax.set_yticks(ticks=[0.0, 0.25*np.sqrt(3.0), 0.5*np.sqrt(3.0)], labels='') # 2次元座標のy軸目盛 ax.set_zlabel('density') # z軸ラベル ax.set_box_aspect(aspect=(1, 1, 1)) # アスペクト比 ax.set_title(label=param_text, loc='left') # パラメータラベル ax.view_init(elev=40, azim=h) # 表示角度 # gif画像を作成 ani = FuncAnimation(fig=fig, func=update, frames=frame_num, interval=100) # gif画像を保存 ani.save('../figure/ternary_3d_turn.gif')
各フレームの作図処理を関数として定義して、FuncAnimation()を使ってアニメーション(gif画像)を作成する。
ax.view_init()で表示角度を設定できる。フレームごとに水平方向の角度の引数azimに指定する値を変更する。
ディリクレ分布のパラメータとして利用する値を指定する。
# ディリクレ分布のパラメータとして利用する値を指定 alpha_0_vals = np.arange(start=1.0, stop=10.1, step=0.1).round(decimals=1) alpha_1_vals = np.arange(start=2.0, stop=11.1, step=0.1).round(decimals=1) alpha_2_vals = np.arange(start=3.0, stop=12.1, step=0.1).round(decimals=1) # フレーム数を設定 frame_num = len(alpha_0_vals) # z軸の最小値と最大値を設定 dens_min = 0.0 dens_max = 27.0 alpha_max_k = np.array([alpha_0_vals.max(), alpha_1_vals.max(), alpha_2_vals.max()]) dens_max = np.ceil( dirichlet.pdf(x=(alpha_max_k-1.0)/(np.sum(alpha_max_k)-3.0), alpha=alpha_max_k) ) # 等高線を引く値を設定 dens_levels = np.linspace(dens_min, dens_max, num=10) print(dens_levels)
[ 0. 3. 6. 9. 12. 15. 18. 21. 24. 27.]
詳しくは前回の記事を参照のこと。
曲面図のアニメーションを作成する。
# 三角座標上の等高線図を作成 fig = plt.figure(figsize=(12, 12), facecolor='white') # 図の設定 ax = fig.add_subplot(projection='3d') # 3D用の設定 fig.suptitle(t='3D Ternary Plot', fontsize=20) # タイトル # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目のパラメータを取得 alpha_k = np.array([alpha_0_vals[i], alpha_1_vals[i], alpha_2_vals[i]]) # ディリクレ分布の確率密度を計算 dens_vals = np.array( [dirichlet.pdf(x=x_k, alpha=alpha_k) if all(x_k != np.nan) else np.nan for x_k in x_points] ) # パラメータラベル用の文字列を作成 param_text = '$' + '\\alpha=('+', '.join([str(val) for val in alpha_k])+')' + ', x=(x_0, x_1, x_2)' + '$' # 三角座標上の曲面図を作成 ax.quiver(grid_x, grid_y, np.zeros_like(grid_x), grid_u, grid_v, np.zeros_like(grid_x), arrow_length_ratio=0.0, ec='gray', linewidth=1.5, linestyle=':') # 三角座標のグリッド線 ax.quiver(axis_x, axis_y, np.zeros_like(axis_x), axis_u, axis_v, np.zeros_like(axis_x), arrow_length_ratio=0.0, ec='black', linestyle='-') # 三角座標の枠線 for val in axis_vals: ax.text(x=0.5*val-0.05, y=0.5*val*np.sqrt(3.0), z=0.0, s=str(np.round(1.0-val, 1)), ha='center', va='center') # 三角座標のx軸目盛 ax.text(x=val, y=0.0-0.05, z=0.0, s=str(np.round(val, 1)), ha='center', va='center') # 三角座標のy軸目盛 ax.text(x=0.5*val+0.5+0.05, y=0.5*(1.0-val)*np.sqrt(3.0), z=0.0, s=str(np.round(1.0-val, 1)), ha='center', va='center') # 三角座標のz軸目盛 ax.text(x=0.25-0.1, y=0.25*np.sqrt(3.0), z=0.0, s='$x_0$', ha='right', va='center', size=25) # 三角座標のx軸ラベル ax.text(x=0.5, y=0.0-0.1, z=0.0-0.1, s='$x_1$', ha='center', va='top', size=25) # 三角座標のy軸ラベル ax.text(x=0.75+0.1, y=0.25*np.sqrt(3.0), z=0.0, s='$x_2$', ha='left', va='center', size=25) # 三角図のz軸ラベル ax.contour(y_0_grid, y_1_grid, dens_vals.reshape(y_shape), vmin=dens_min, vmax=dens_max, levels=dens_levels, offset=0.0) # 確率密度の等高線 ax.plot_surface(y_0_grid, y_1_grid, dens_vals.reshape(y_shape), cmap='viridis', alpha=0.8) # 確率密度の曲面 ax.set_xticks(ticks=[0.0, 0.5, 1.0], labels='') # 2次元座標のx軸目盛 ax.set_yticks(ticks=[0.0, 0.25*np.sqrt(3.0), 0.5*np.sqrt(3.0)], labels='') # 2次元座標のy軸目盛 ax.set_zlabel('density') # z軸ラベル ax.set_zlim(bottom=dens_min, top=dens_max) # z軸の表示範囲 ax.set_box_aspect(aspect=(1, 1, 1)) # アスペクト比 ax.set_title(label=param_text, loc='left') # パラメータラベル # gif画像を作成 ani = FuncAnimation(fig=fig, func=update, frames=frame_num, interval=100) # gif画像を保存 ani.save('../figure/ternary_3d_dens.gif')
各フレームの確率密度の計算と作図処理を関数として定義して、アニメーションを作成する。
以上で、3次元の三角図を作成できた。
おわりに
以上、ディリクレ分布の記事を書くために記事が3つ生えました。
2022年10月18日はJuice=Juiceの入江里咲さんの17歳のお誕生日です！
かわいい///
