Zenn
🟠

多変量ラプラス分布からのデータ生成

2023/11/17に公開
Python
math
#
統計
#
統計学
#
データサイエンス
tech

はじめに

論文に載っているシミュレーションのデータを生成しようと思ったら、多変量ラプラス分布がでてきて、scipyで実装しようと思ったら、まさかのなくて自分で関数作らないといけなくなったので、もし同じ状況になった方いたら参考にしてください。

内容は詳しくはこの論文を見てください。(参考文献にリンクあります)
The Laplace Distribution and Generalizations

多変量ラプラス分布

密度関数

平均0 、共分散行列が\Sigmaのラプラス分布の密度関数は次のようにかける。

g(\bm{y}) = \frac{2}{(2\pi)^{d/2}|\Sigma|^{1/2}}\left( \frac{\bm{y}^T\Sigma^{-1}\bm{y}}{2} \right)^{\nu/2}K_{\nu}\left( \sqrt{2\bm{y}'\Sigma^{-1}\bm{y}} \right),

ここで v = (2 − d)/2

K_\nu(*) は次の(1)または(2)である。

\begin{align} K_{\lambda}(u) &= \frac{1}{2} \left(\frac{u}{2}\right)^{\lambda} \int_{0}^{\infty} t^{-\lambda-1} \exp\left(-t - \frac{u^2}{4t}\right) dt, \quad u > 0. \\ K_{\lambda}(u) &= \frac{(u/2)^{\lambda}\Gamma(1/2)}{\Gamma(\lambda+1/2)} \int_{1}^{\infty} e^{-ut}(t^2-1)^{\lambda-1/2} dt, \quad \lambda \geq -1/2. \end{align}

この定義見て、ガンマ関数使って式変形すれば(1)が(2)になるのかなーって思ったけど、使い方としては、d=2 , 3のとき(2)でその他は(1)使う感じかな。

近似

実際、上記の密度関数を使ってpythonで実装するのはかなりめんどくさいですが、とても便利な近似があるので実装はこちらを使います。

X \sim N_d(0 , \Sigma)とする。つまりXは平均0,共分散行列\Sigmaをもつd次元正規分布から発生。

W \sim Exp(1) とする。つまりWは平均1,分散1の指数分布から発生。Xとは独立。

Y=\sqrt W Xと変換する。

するとこのYが平均0,共分散行列\Sigmaのラプラス分布に近似できる。

E(Y)=E(\sqrt W)E(X)=0

Var(Y)=E(YY^T)=E(W)E(XX^T) =\Sigma

Pythonで実装

まずはライブラリをインストール

import numpy as np
from scipy.stats import multivariate_normal, expon
import matplotlib.pyplot as plt

多変量ラプラス分布を定義する。

def multivariate_laplace(mean, covariance, size):
    # d次元正規分布からサンプルを生成
    x = multivariate_normal.rvs(mean=mean, cov=covariance, size=size)

    # 指数分布からサンプルを生成し、正規分布のサンプルに掛け合わせる
	  W = expon.rvs(scale=1, size=size)
    return \sqrt(W[:, np.newaxis]) * x

データを生成する。

mean = [0, 0]  # 平均
covariance = [[1, 0], [0, 1]]  # 共分散行列

# サンプルの生成
samples = multivariate_laplace(mean, covariance, size=1000)

ここまででデータの生成はひとまず終わりです。

最後にラプラス分布と正規分布をプロットして比較してみる。

# 平均と共分散行列を設定
mean = [0, 0]  # 平均
covariance = [[1, 0], [0, 1]]  # 共分散行列

# サンプルの生成
samples_laplace = multivariate_laplace(mean, covariance, size=1000)
samples_normal = multivariate_normal.rvs(mean, covariance, size=1000)

# プロットの作成
fig, axs = plt.subplots(1, 2, figsize=(12, 6), sharex=True, sharey=True)

# ラプラス分布から生成したサンプルのプロット
axs[0].scatter(samples_laplace[:, 0], samples_laplace[:, 1], alpha=0.5)
axs[0].set_title('Multivariate Laplace Distribution')
axs[0].set_xlim([-5, 5])  # x軸の範囲を設定
axs[0].set_ylim([-5, 5])  # y軸の範囲を設定

# 標準正規分布から生成したサンプルのプロット
axs[1].scatter(samples_normal[:, 0], samples_normal[:, 1], alpha=0.5, color='orange')
axs[1].set_title('Multivariate Normal Distribution')
axs[1].set_xlim([-5, 5])  # x軸の範囲を設定
axs[1].set_ylim([-5, 5])  # y軸の範囲を設定

plt.tight_layout()
plt.show()

ラプラス分布の方が中心に集まりつつも裾に広がっている印象です。

まとめ

詳しいことは参考文献に載っている論文を見てください。平均\muのラプラス分布や非対称のラプラス分布も扱っています。

参考文献

Discussion