Help us understand the problem. What is going on with this article?

PRMLの図作成 演習問題1.4 確率密度関数の非線形変換

More than 3 years have passed since last update.

やりたいこと

確率変数$x$と$y$が、$x = g(y)$の関係になっているとします。また、$x$に関する確率密度関数$p_x(x)$は分かっているとします。
この時、$y$に関する確率密度関数$p_y(y)$の形はどうなっているかを考えます。

説明

単純に変換すると、$p_x(g(y))$となりますが、これは$x$に関する確率密度関数であるため、$p_y(y) \neq p_x(g(y))$です。

$p_x(x)$と$p_y(y)$の関係は、任意の範囲$x_1\sim x_2$に関して、

\int_{x_1}^{x_2} p_x(x) \mathrm{d}x = \int_{g^{-1}(x_1)}^{g^{-1}(x_2)} p_y(y) \mathrm{d}y

となるはずです。
ただし、$g^{-1}(x)$は、$g(x)$の逆関数です。

積分の変数変換の公式を用いると、

\begin{align}
\int_{x_1}^{x_2} p_x(x) \mathrm{d}x&=\int_{x_1}^{x_2} p_x(g(y)) \mathrm{d}x\\
&=\int_{g^{-1}(x_1)}^{g^{-1}(x_2)} p_x(g(y))\left|
\frac{\partial g(y)}{\partial y}\right| \mathrm{d}y
\end{align}

よって、

p_y(y) = p_x(g(y))\left|
\frac{\partial g(y)}{\partial y}\right|

となります。

実装

p_x(x) = \mathcal{N}(x\mid\mu,\sigma^2)\\
x = g(y) = \ln(y)-\ln(1-y)+5\\
y = g^{-1}(x) = \frac{1}{1+\exp(-x+5)}

の場合を考えます。

\frac{\partial g(y)}{\partial y} = \frac{1}{y(1-y)}

より、

p_y(y) = \mathcal{N}(g(y)\mid\mu,\sigma^2)\frac{1}{y(1-y)} 

これが、$p_x(x)$からのサンプルを$y=g^{-1}(x)$で変換したデータのヒストグラムと一致するかを確認します。

コード

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

#ガウス分布の密度関数
def gaussianDist(sig,mu,x):
    y=np.exp(-(x-mu)**2/(2*sig**2))/(np.sqrt(2*np.pi)*sig)
    return y

#確率変数の変換関数
def g(y):
    x=np.log(y)-np.log(1-y)+5
    return x

#確率変数の変換関数の逆関数
def invg(x):
    y=1/(1+np.exp(-x+5))
    return y

#ガウス分布px(x)の平均、分散
sig=1.0
mu=6

#ヒストグラムのサンプル数
N = 50000 

plt.xlim([0,10])
plt.ylim([0,1])

####
x = np.linspace(0,10,100)

#確率変数の変換関数をプロット
y=invg(x)
plt.plot(x,y,'b')

#px(x)のプロット
y = gaussianDist(sig,mu,x)
plt.plot(x,y,'r')

#px(x)からのサンプルを元にヒストグラムをプロット
x_sample = mu + sig * np.random.randn(N)
plt.hist(x_sample,bins=20,normed=True,color='lavender')


####
y=np.linspace(0.01,0.99,100)

##py(y)のプロット
x=gaussianDist(sig,mu,g(y))/(y*(1-y))
plt.plot(x,y,'m')

#px(x)からのサンプルをg^-1(x)で変換したデータのヒストグラムをプロット
y_sample = invg(mu + sig * np.random.randn(N))
plt.hist(y_sample,bins=20,normed=True,orientation="horizontal",color='lavender')

#px(g(y))のように単純に変換した関数をプロット
x = gaussianDist(sig,mu,g(y))
plt.plot(x/(x.sum()*0.01) ,y,'lime')

####
#平均muとg^-1(mu)との関係をプロット
plt.plot([mu, mu], [0, invg(mu)], 'k--')
plt.plot([0, mu], [invg(mu), invg(mu)], 'k--')

実行結果

test.png

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away