使用するモジュール

matplotlib, numpy

コード書くとき参考にした文献

確率楕円のみ生成するコード

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

###################### 設定 ######################
n_sigma = 3 # nシグマ(標準偏差)分を覆う．3の場合は3σ確率楕円を生成．

###################### データの準備 ######################
data = np.random.randn(1000,2) # 100x2行列で，標準正規分布の乱数を生成　

###################### 確率楕円の生成 ######################
x = data[:,0]; y = data[:,1]

cov = np.cov(x, y) # 共分散
eigenvalues, eigenvectors = np.linalg.eig(cov) # 固有値と固有ベクトル

w, h = 2 * n_sigma * np.sqrt(eigenvalues) # 長軸と短軸の長さ
theta_rad = np.arctan2(*eigenvectors[:,0][::-1]) # 楕円の傾き，radian
theta = np.degrees(theta_rad) # 楕円の傾き，degree

###################### グラフの生成 ######################
ax = plt.subplot(111)
ax.set_aspect('equal')
ax.scatter(x, y)

ell = Ellipse(xy=(np.mean(x), np.mean(y)),width=w, height=h, angle=theta, color='black')
ell.set_facecolor('none')
ax.add_artist(ell)

plt.show()

結果

序に確率楕円の座標を取得し，外接する長方形も表示するコード

import numpy as np
import matplotlib.pyplot as plt

###################### 設定 ######################
n_sigma = 3 # nシグマ(標準偏差)分を覆う．3の場合は3σ確率楕円を生成．

###################### データの準備 ######################
data = np.random.randn(1000,2) # 100x2行列で，標準正規分布の乱数を生成　

###################### 確率楕円の生成 ######################
x = data[:,0]; y = data[:,1]

cov = np.cov(x, y) # 共分散
eigenvalues, eigenvectors = np.linalg.eig(cov) # 固有値と固有ベクトル

w, h = 2 * n_sigma * np.sqrt(eigenvalues) # 長軸と短軸の長さ
theta_rad = np.arctan2(*eigenvectors[:,0][::-1]) # 楕円の傾き，radian
theta = np.degrees(theta_rad) # 楕円の傾き，degree

phi = np.linspace(0, 2*np.pi, 1000) # 極座標，角度の列
center_x = np.mean(x); center_y = np.mean(y) # 楕円の中心

rotate = np.array([[np.cos(theta_rad), -np.sin(theta_rad)],[np.sin(theta_rad), np.cos(theta_rad)]]) # 回転する行列
ellipse = np.array([w/2*np.cos(phi), h/2*np.sin(phi)]) # 回転，移動させる前の楕円
ellipse = (np.matmul(rotate, ellipse)+np.array([[center_x],[center_y]])).T # 回転，移動後の楕円
square = np.array([[w/2, -w/2, -w/2, w/2, w/2],[h/2, h/2, -h/2, -h/2, h/2]]) # 回転，移動させる前の長方形
square = (np.matmul(rotate, square)+np.array([[center_x],[center_y]])).T # 回転，移動後の長方形

###################### グラフの生成 ######################
ax = plt.subplot(111)
ax.set_aspect('equal')
ax.scatter(x, y)
ax.plot(ellipse[:,0], ellipse[:,1], '--',lw=1)
ax.plot(square[:,0], square[:,1], '--',lw=1)

plt.show()

結果