[PYTHON] PRML Diagram Drawing Exercise 1.4 Nonlinear Transformation of Probability Density Function

Thing you want to do

Suppose the random variables $ x $ and $ y $ have a relationship of $ x = g (y) $. We also know the probability density function $ p_x (x) $ for $ x $. At this time, consider what the shape of the probability density function $ p_y (y) $ for $ y $ is.

Description

A simple conversion would be $ p_x (g (y)) $, which is a probability density function for $ x $, so it is $ p_y (y) \ neq p_x (g (y)) $.

The relationship between $ p_x (x) $ and $ p_y (y) $ is for any range $ 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

Should be. However, $ g ^ {-1} (x) $ is the inverse function of $ g (x) $.

Using the change of variables formula of the integral,

\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}

Therefore,

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

It will be.

Implementation

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)}

Consider the case of.

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

Than,

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

Check if this matches the histogram of the data from $ p_x (x) $ converted with $ y = g ^ {-1} (x) $.

code

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

import numpy as np
import matplotlib.pyplot as plt

#Gaussian distribution density function
def gaussianDist(sig,mu,x):
    y=np.exp(-(x-mu)**2/(2*sig**2))/(np.sqrt(2*np.pi)*sig)
    return y

#Random variable conversion function
def g(y):
    x=np.log(y)-np.log(1-y)+5
    return x
    
#Inverse function of random variable conversion function
def invg(x):
    y=1/(1+np.exp(-x+5))
    return y

#Gaussian distribution px(x)Mean, variance
sig=1.0
mu=6

#Histogram sample size
N = 50000 

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

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

#Plot the random variable transformation function
y=invg(x)
plt.plot(x,y,'b')

#px(x)Plot
y = gaussianDist(sig,mu,x)
plt.plot(x,y,'r')

#px(x)Plot the histogram based on the sample from
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)Plot
x=gaussianDist(sig,mu,g(y))/(y*(1-y))
plt.plot(x,y,'m')

#px(x)Sample from g^-1(x)Plot the histogram of the data converted in
y_sample = invg(mu + sig * np.random.randn(N))
plt.hist(y_sample,bins=20,normed=True,orientation="horizontal",color='lavender')

#px(g(y))Plot the converted function simply like
x = gaussianDist(sig,mu,g(y))
plt.plot(x/(x.sum()*0.01) ,y,'lime')

####
#Mean mu and g^-1(mu)Plot the relationship with
plt.plot([mu, mu], [0, invg(mu)], 'k--')
plt.plot([0, mu], [invg(mu), invg(mu)], 'k--')

Execution result