2変量正規分布の3D描画
Pythonスクリプト
2変量正規分布は、各変量の分散が1、相関係数がであるとき、次式で表される。
式(1)において、の場合の描画例が図1である。
図1
スクリプトは以下に示すとおりである。
============================================
from mpl_toolkits.mplot3d import
Axes3D
import matplotlib.pyplot as plt
import numpy as np
import math
class BiVarNorm:
"""
Bivariate normal distribution
"""
def __init__(self, rho):
self.rho = rho
def value(self, x, y):
v =
(x**2 + y**2 - 2 * self.rho * x * y) / (2.0 * (1.0 - self.rho**2))
v =
np.exp(-v) / ((2.0 * math.pi * (1.0 - self.rho**2)) ** 0.5)
return v
while True:
rho = float(input('rho = '))
if math.fabs(rho) >= 1.0:
break
fig = plt.figure()
ax = Axes3D(fig)
Xg = np.arange(-4.0, 4.0,
0.1)
Yg = np.arange(-4.0, 4.0,
0.1)
X, Y = np.meshgrid(Xg, Yg)
bi_norm = BiVarNorm(rho)
ax.plot_surface(X, Y,
bi_norm.value(X, Y), rstride = 1, cstride = 1, cmap = plt.cm.coolwarm)
plt.show()
print('\n終了です。\n')
============================================
次の関数valueにおいて
def value(self, x, y):
v =
(x**2 + y**2 - 2 * self.rho * x * y) / (2.0 * (1.0 - self.rho**2))
v =
np.exp(-v) / ((2.0 * math.pi * (1.0 - self.rho**2)) ** 0.5)
return v
指数関数np.expが用いられているが、これをmath.expに変えると「配列に対応していない」旨のエラーが表示される。配列(リスト)を引数とする関数を用いずに関数ax.plot_surfaceの第3引数用の配列を用意するスクリプト例は別のウェブサイトに用意している。
上のスクリプトを実行すると「rho = 」との値の入力が求められる。
図2
絶対値が1より小さい適当な値を入力すると、その値を相関係数とする式(1)のグラフが描かれる(図1)。描画画像上をマウスでドラッグすると、画像を回転することができる(図3)。
図3
表示されている画像のフォームの右上のX印をクリックしてフォームを閉じると、次の相関係数の値の設定が求められる(図4)。
図4
絶対値が1より小さい値を設定すると、その値を相関係数とするグラフが描かれる。1以上の値を設定すれば、プログラムの実行終了となる。
プログラム(スクリプト)は圧縮ファイルPrg.zipとして用意した。ファイル名Prg.zipをクリックしてダウンロードすることができる。