他人に PCA（主成分分析）などを説明する際に、いつも「楕円形に分布する乱数を作りたいなぁ」と思う。
あるあるですね。
でもその度に毎回
- 各軸がガウス分布に従う２次元の変数だから、numpy.random の multivariate_normal か
- cov ってなんだっけ・・・、あ、共分散行列か
- 作りたい形の共分散行列ってどうやって作るんだよ・・・
- （適当な値の cov をつっこむ）"covariance is not symmetric positive-semidefinite." いやまぁそうだろうけどね
となって、途方に暮れる。
なので、ちゃんと作って見た。
関数
from numpy.random import Generator, default_rng
from numpy.typing import NDArray
def generate_random_ellipse(
corrcoef: float,
degree: float,
size: int,
seed: int = 12,
) -> NDArray:
"""楕円形に分布する２次元の乱数を生成する。
Args:
corrcoef (float):
生成する乱数の、角度を回転させる前の相関係数。
楕円の平べったさになる（回転した後の相関係数はこの値とは異なる）。
0 だと真ん丸で、1 だと平べったい。
degree (float):
楕円の長径の x 軸との角度。
size (float):
乱数の個数。
seed (int):
乱数シード。
Returns:
NDArray: (size, 2) の行列
"""
rng: Generator = default_rng(seed)
x: NDArray = rng.normal(loc=0, scale=1, size=size)
y: NDArray = rng.normal(loc=0, scale=1, size=size)
z: NDArray = corrcoef * x + (1 - (corrcoef**2)) ** 0.5 * y
d: NDArray = np.concatenate([x.reshape(1, -1), z.reshape(1, -1)], axis=0)
rad: float = np.deg2rad(degree - 45)
m: NDArray = np.array(
[[np.cos(rad), -np.sin(rad)], [np.sin(rad), np.cos(rad)]]
)
return np.dot(m, d).T
相関のある２つの疑似乱数の生成は、こちらを参照させていただきました。
上記で楕円形の乱数を発生させて、それを回転させているだけ。
使ってみる
例えばこんな感じで生成します。
x: NDArray = generate_random_ellipse(corrcoef=0.8, degree=70, size=1000)
するとこんなのが出来ます。
もういつでも楕円形に分布する乱数が作れますね！
ちなみに
作った乱数の相関係数は、もちろん指定した値とは変わってしまいます。
>>> np.corrcoef(x, rowvar=False)
array([[1. , 0.63088707],
[0.63088707, 1. ]])
そしてこの乱数の共分散行列は以下。
>>> np.cov(x, rowvar=False)
array([[0.38510963, 0.49677686],
[0.49677686, 1.61003414]])
これを numpy.random の multivariate_normal にあげれば、同様の乱数が作れます。
rng: Generator = default_rng(12)
x2: NDArray = rng.multivariate_normal(
mean=np.array([0, 0]),
cov=np.array([[0.38510963, 0.49677686], [0.49677686, 1.61003414]]),
size=1000,
)
