はじめに
PythonのMatplotlibライブラリを使ってグラフを動かそうシリーズです。
この記事では、全ての方向に移動する3次元のランダムウォークのアニメーションを作成します。
【前の内容】
【他の内容】
【目次】
3次元ランダムウォークの作図：全方向移動
空間上の全ての方向(360°×360°)に移動するランダムウォーク(random walk)のアニメーションを作成します。
利用するライブラリを読み込みます。
# 利用するモジュール import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation import numpy as np
全方向移動の確認
まずは、空間上を移動する処理を確認します。
球面の座標計算やラジアンについては「球面の作図【Matplotlib】 - からっぽのしょこ」を参照してください。
・作図コード(クリックで展開)
単位円の円周上の点の座標を計算します。
# フレーム数を指定 frame_num = 101 # 点用のラジアンを作成 point_t_vec = np.linspace(start=0, stop=2*np.pi, num=frame_num+1)[:frame_num+1] # 1周期分 point_u_vec = np.linspace(start=0, stop=2*np.pi, num=frame_num+1)[:frame_num+1] # 1周期分 print(point_t_vec[:5].round(2)) print(point_u_vec[:5].round(2)) # 単位球面上の点の座標を計算 point_x_vec = np.sin(point_t_vec) * np.cos(point_u_vec) point_y_vec = np.sin(point_t_vec) * np.sin(point_u_vec) point_z_vec = np.cos(point_t_vec) print(point_x_vec[:5].round(2)) print(point_y_vec[:5].round(2)) print(point_z_vec[:5].round(2))
[0. 0.06 0.12 0.19 0.25]
[0. 0.06 0.12 0.19 0.25]
[0. 0.06 0.12 0.18 0.24]
[0. 0. 0.02 0.03 0.06]
[1. 1. 0.99 0.98 0.97]
フレーム数 frame_num を指定して、frame_num 個の点の座標を作成します。
から
の値
を作成して、x軸の値を
、y軸の値を
、z軸の値を
で計算します。
は円周率で、
np.pi で扱えます。
ノルムが1の目盛線(単位球の円周)の描画用の値を作成します。
# 球面用のラジアンを作成 sphere_t_vec = np.linspace(start=0.0, stop=2.0*np.pi, num=61) sphere_u_vec = np.linspace(start=0.0, stop=2.0*np.pi, num=61) sphere_t_arr, sphere_u_arr = np.meshgrid(sphere_t_vec, sphere_u_vec) # 格子点を作成 print(sphere_t_arr[:3, :5].round(2)) print(sphere_u_arr[:3, :5].round(2)) # 球面の座標を計算 sphere_x_arr = np.sin(sphere_t_arr) * np.cos(sphere_u_arr) sphere_y_arr = np.sin(sphere_t_arr) * np.sin(sphere_u_arr) sphere_z_arr = np.cos(sphere_t_arr) print(sphere_x_arr[:3, :5].round(2)) print(sphere_y_arr[:3, :5].round(2)) print(sphere_z_arr[:3, :5].round(2))
[[0. 0.1 0.21 0.31 0.42]
[0. 0.1 0.21 0.31 0.42]
[0. 0.1 0.21 0.31 0.42]]
[[0. 0. 0. 0. 0. ]
[0.1 0.1 0.1 0.1 0.1 ]
[0.21 0.21 0.21 0.21 0.21]]
[[0. 0.1 0.21 0.31 0.41]
[0. 0.1 0.21 0.31 0.4 ]
[0. 0.1 0.2 0.3 0.4 ]]
[[0. 0. 0. 0. 0. ]
[0. 0.01 0.02 0.03 0.04]
[0. 0.02 0.04 0.06 0.08]]
[[1. 0.99 0.98 0.95 0.91]
[1. 0.99 0.98 0.95 0.91]
[1. 0.99 0.98 0.95 0.91]]
2つのラジアンを格子状に複製して、球面の座標を作成します。
作成した点を描画します。
# グラフサイズを設定 axis_size = 1.5 # グラフオブジェクトを初期化 fig, ax = plt.subplots(figsize=(9, 9), facecolor='white', subplot_kw={'projection': '3d'}) fig.suptitle('unit spherical surface', fontsize=20) # 矢サイズを指定 alr = 0.1 # 作図処理を定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # 点の座標を取得 t = point_t_vec[i] u = point_u_vec[i] x = point_x_vec[i] y = point_y_vec[i] z = point_z_vec[i] # ラベル用の文字列を作成 coord_label = '$(t, u) = ({:.2f}\pi, {:.2f}\pi), (x, y, z) = ({:.2f}, {:.2f}, {:.2f})$'.format( t/np.pi, u/np.pi, x, y, z ) # 単位球面状の点を作図 ax.quiver([-axis_size, 0, 0], [0, -axis_size, 0], [0, 0, -axis_size], [2*axis_size, 0, 0], [0, 2*axis_size, 0], [0, 0, 2*axis_size], arrow_length_ratio=alr/axis_size, color='black', linewidth=1) # 軸線 ax.plot_wireframe(sphere_x_arr, sphere_y_arr, sphere_z_arr, color='gray', alpha=0.1) # 球面 ax.plot([0, x], [0, y], [0, z], linewidth=2) # 動径 ax.scatter(x, y, z, s=100) # 球面上の点 ax.plot(point_x_vec[:i+1], point_y_vec[:i+1], point_z_vec[:i+1], linewidth=1) # 軌跡 ax.set_xlim(xmin=-axis_size, xmax=axis_size) ax.set_ylim(ymin=-axis_size, ymax=axis_size) ax.set_zlim(zmin=-axis_size, zmax=axis_size) ax.set_xlabel('$x = \sin t\ \cos u$') ax.set_ylabel('$y = \sin t\ \sin u$') ax.set_zlabel('$z = \cos t$') ax.set_title(coord_label, loc='left') ax.set_box_aspect((1, 1, 1)) ax.view_init(elev=30, azim=i) # 表示アングル # 動画を作成 ani = FuncAnimation(fig=fig, func=update, frames=frame_num, interval=100) # 動画を書出 ani.save(filename='UnitSphere.gif', dpi=100)
2つのラジアン(弧度法の角度) 、
を用いて、単位球(半径が1)の球面上の点の座標
は、次の式で得られます。
は円周率です。
つまり、 から
の乱数を2つ生成して座標を計算すると、現在の位置(球の中心)からランダムな方向に距離(ノルム)1の移動を表せます。
1サンプル
次は、1つのサンプルを描画します。
試行回数を指定して、移動量の計算用の値をランダムに生成します。
# 試行回数を指定 max_iter = 100 # 乱数を生成 random_t_vec = 2*np.pi * np.random.uniform(low=0, high=1, size=max_iter) # ラジアン random_u_vec = 2*np.pi * np.random.uniform(low=0, high=1, size=max_iter) # ラジアン print(random_t_vec[:5].round(2)) print(random_t_vec.shape) print(random_u_vec[:5].round(2)) print(random_u_vec.shape)
[5.2 3.18 3.68 5.76 3.61]
(100,)
[3.65 4.55 1.7 5.37 1.22]
(100,)
から
の一様乱数を
np.random.uniform()で生成します。 は円周率で、
np.pi で扱えます。最小値の引数 low に 0、最大値の引数 high に 1 を指定して、2*pi を掛けます。または、high 引数に 2*pi を指定します。
この例では原理的な処理が分かりやすいように、0 から 1 の乱数に 2*np.pi を掛けることで から
の乱数を作成しています。
移動量に応じた座標を計算します。
# 試行ごとに集計 x_vec = np.cumsum( np.hstack([0, np.sin(random_t_vec)*np.cos(random_u_vec)]) # 初期値を追加 ) y_vec = np.cumsum( np.hstack([0, np.sin(random_t_vec)*np.sin(random_u_vec)]) # 初期値を追加 ) z_vec = np.cumsum( np.hstack([0, np.cos(random_t_vec)]) # 初期値を追加 ) print(x_vec[:5].round(2)) print(x_vec.shape) print(y_vec[:5].round(2)) print(y_vec.shape) print(z_vec[:5].round(2)) print(z_vec.shape)
[0. 0.77 0.78 0.85 0.54]
(101,)
[ 0. 0.43 0.47 -0.04 0.36]
(101,)
[ 0. 0.47 -0.53 -1.39 -0.52]
(101,)
乱数を用いて各座標を np.cos() と np.sin() を使って計算します。
初期値(スタート地点・0回目の結果)を 0 として追加して、各試行までの和(累積和)を np.cumsum() で計算します。
点の推移のグラフを作成します。
# 最終値を取得 x = x_vec[max_iter] y = y_vec[max_iter] z = z_vec[max_iter] # 余白サイズを指定 margin_rate = 0.1 # グラフサイズを設定 axis_size = np.max([np.abs(x_vec), np.abs(y_vec), np.abs(z_vec)]) axis_size = np.ceil(axis_size * (1+margin_rate)) # 余白を追加 # 矢サイズを指定 alr = 0.5 # 3Dランダムウォークを作図 fig, ax = plt.subplots(figsize=(12, 12), dpi=250, facecolor='white', subplot_kw={'projection': '3d'}) ax.quiver([-axis_size, 0, 0], [0, -axis_size, 0], [0, 0, -axis_size], [2*axis_size, 0, 0], [0, 2*axis_size, 0], [0, 0, 2*axis_size], arrow_length_ratio=alr/axis_size, color='black', linewidth=1) # x・y・z軸線 ax.plot(x_vec, y_vec, zs=-axis_size, zdir='z', color='C0', linewidth=1, linestyle='dashed') # xy面の軌跡 ax.plot(y_vec, z_vec, zs=-axis_size, zdir='x', color='C0', linewidth=1, linestyle='dashed') # yz面の軌跡 ax.plot(x_vec, z_vec, zs=axis_size, zdir='y', color='C0', linewidth=1, linestyle='dashed') # xz面の軌跡 ax.scatter([x, -axis_size, x], [y, y, axis_size], [-axis_size, z, z], s=100, fc='none', ec='blue') # 面ごとの最終地点 ax.quiver(x, y, z, [0, -axis_size-x, 0], [0, 0, axis_size-y], [-axis_size-z, 0, 0], arrow_length_ratio=0, color='gray', linestyle='dotted') # 点の目盛線 ax.plot(x_vec, y_vec, z_vec) # 軌跡 ax.scatter(x, y, z, s=100, color='blue') # 最終地点 ax.margins(x=0.0, y=0.0, z=0.0) # (機能しない?のでlimで謎小細工) ax.set_xlim(xmin=-axis_size/1.05, xmax=axis_size/1.05) ax.set_ylim(ymin=-axis_size/1.05, ymax=axis_size/1.05) ax.set_zlim(zmin=-axis_size/1.05, zmax=axis_size/1.05) ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') ax.set_title('iteration: {}, $(x, y, z) = ({:.2f}, {:.2f}, {:.2f})$'.format(max_iter, x, y, z), loc='left') fig.suptitle('Random Walk', fontsize=20) ax.set_box_aspect((1, 1, 1)) plt.show()
各軸線の最小値・最大値(各方向のグラフサイズ)は Axes.set_*lim() に指定した値に余白(マージン)が追加されます。よって、axis_size の位置にプロットすると側面や底面からズレます。マージン率(マージンを含めた全体のサイズは (最大値 - 最小値) * 1.05 らしい)を0にできればいいのですがやり方が分からなかったので、回避策としてごにょごにょしています(まっとうなやり方を教えてください)。コードの分かりやすさを優先するのであれば 1/1.05 を外してください。
点の軌道(実線)と軸ごとの軌道(破線)が近いと分かりづらいので margin_rate で調整しています。
点の推移のアニメーションを作成します。
# 余白サイズを指定 margin_rate = 0.1 # グラフサイズを設定 axis_size = np.max([np.abs(x_vec), np.abs(y_vec), np.abs(z_vec)]) axis_size = np.ceil(axis_size * (1+margin_rate)) # 余白を追加 # グラフオブジェクトを初期化 fig, ax = plt.subplots(figsize=(12, 12), facecolor='white', subplot_kw={'projection': '3d'}) fig.suptitle('Random Walk', fontsize=20) # 矢サイズを指定 alr = 0.5 # 作図処理を定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # 現在値を取得 x = x_vec[i] y = y_vec[i] z = z_vec[i] # 3Dランダムウォークを作図 ax.quiver([-axis_size, 0, 0], [0, -axis_size, 0], [0, 0, -axis_size], [2*axis_size, 0, 0], [0, 2*axis_size, 0], [0, 0, 2*axis_size], arrow_length_ratio=alr/axis_size, color='black', linewidth=1) # x・y・z軸線 ax.plot(x_vec[:i+1], y_vec[:i+1], zs=-axis_size, zdir='z', color='C0', linewidth=1, linestyle='dashed') # xy面の軌跡 ax.plot(y_vec[:i+1], z_vec[:i+1], zs=-axis_size, zdir='x', color='C0', linewidth=1, linestyle='dashed') # yz面の軌跡 ax.plot(x_vec[:i+1], z_vec[:i+1], zs=axis_size, zdir='y', color='C0', linewidth=1, linestyle='dashed') # xz面の軌跡 ax.scatter([x, -axis_size, x], [y, y, axis_size], [-axis_size, z, z], s=100, fc='none', ec='blue') # 面ごとの現在地点 ax.quiver(x, y, z, [0, -axis_size-x, 0], [0, 0, axis_size-y], [-axis_size-z, 0, 0], arrow_length_ratio=0, color='gray', linestyle='dotted') # 点の目盛線 ax.plot(x_vec[:i+1], y_vec[:i+1], z_vec[:i+1]) # 軌跡 ax.scatter(x, y, z, s=100, color='blue') # 現在地点 ax.margins(x=0.0, y=0.0, z=0.0) # (機能しない?のでlimで謎小細工) ax.set_xlim(xmin=-axis_size/1.05, xmax=axis_size/1.05) ax.set_ylim(ymin=-axis_size/1.05, ymax=axis_size/1.05) ax.set_zlim(zmin=-axis_size/1.05, zmax=axis_size/1.05) ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') ax.set_title('iteration: {}, $(x, y, z) = ({:.2f}, {:.2f}, {:.2f})$'.format(i, x, y, z), loc='left') ax.set_box_aspect((1, 1, 1)) # 動画を作成 ani = FuncAnimation(fig=fig, func=update, frames=max_iter+1, interval=100) # 動画を書出 ani.save('3d_360_xyz_1.gif', dpi=100)
2つの軸の変化を組み合わせごとに側面に描画しています。
2軸のグラフは、2次元ランダムウォークのように推移していますが、2つの軸に注目すると距離(ノルム)が1以下で変化するので、異なる図です。
複数サンプル
続いて、複数のサンプルを並べて描画します。
作図コードについては「RandomWalk_3d_360.py at anemptyarchive/matplotlib · GitHub」を参照してください。
・座標軸のみ
等確率で移動する場合、初期値が期待値になりますが分散を持ち(散らばり)ます。
この記事では、全ての方向に移動する3次元ランダムウォークを扱いました。
おわりに
この記事でPython版ランダムウォークシリーズ完結です。想定外に良い感じのグラフを作れて大満足です。
普段は目的なくノリで内容を決めてインプット・アウトプットしているのですが、何のきっかけだったかで書いた球面の作図の内容がここで使えるとは、なんでもやってみるもんですね。
今後も気ままにやっていきます。
【次の内容】
今回の内容の次って何ですか？やはりMCMCか。
