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

test.png

Recommended Posts

PRML Diagram Drawing Exercise 1.4 Nonlinear Transformation of Probability Density Function
Defeat the probability density function of the normal distribution
Precautions when drawing the probability density function and the histogram on top of each other in matplotlib
PRML diagram drawing Figure 1.4 Polynomial curve fitting