[PYTHON] Exercice de dessin de diagramme PRML 1.4 Transformation non linéaire de la fonction de densité de probabilité

Chose que tu veux faire

Supposons que les variables de probabilité $ x $ et $ y $ aient une relation de $ x = g (y) $. Supposons également que vous connaissez la fonction de densité de probabilité $ p_x (x) $ pour $ x $. À ce stade, considérez la forme de la fonction de densité de probabilité $ p_y (y) $ pour $ y $.

La description

Une conversion simple serait $ p_x (g (y)) $, qui est une fonction de densité de probabilité pour $ x $, donc c'est $ p_y (y) \ neq p_x (g (y)) $.

La relation entre $ p_x (x) $ et $ p_y (y) $ est pour toute plage $ 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

Devrait être. Cependant, $ g ^ {-1} (x) $ est la fonction inverse de $ g (x) $.

En utilisant la formule de transformation de variable de l'intégrale,

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

Donc,

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

Ce sera.

la mise en oeuvre

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

Prenons le cas de.

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

Que,

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

Vérifiez si cela correspond à l'histogramme des données de $ p_x (x) $ converties par $ y = g ^ {-1} (x) $.

code

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

import numpy as np
import matplotlib.pyplot as plt

#Fonction de densité de distribution gaussienne
def gaussianDist(sig,mu,x):
    y=np.exp(-(x-mu)**2/(2*sig**2))/(np.sqrt(2*np.pi)*sig)
    return y

#Fonction de conversion de variable probabiliste
def g(y):
    x=np.log(y)-np.log(1-y)+5
    return x
    
#Fonction inverse de la fonction de conversion de la variable stochastique
def invg(x):
    y=1/(1+np.exp(-x+5))
    return y

#Distribution gaussienne px(x)Moyenne, dispersion
sig=1.0
mu=6

#Nombre d'échantillons d'histogramme
N = 50000 

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

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

#Tracer la fonction de conversion des variables stochastiques
y=invg(x)
plt.plot(x,y,'b')

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

#px(x)Tracez l'histogramme en fonction de l'échantillon de
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)Terrain
x=gaussianDist(sig,mu,g(y))/(y*(1-y))
plt.plot(x,y,'m')

#px(x)Échantillon de g^-1(x)Tracez l'histogramme des données converties en
y_sample = invg(mu + sig * np.random.randn(N))
plt.hist(y_sample,bins=20,normed=True,orientation="horizontal",color='lavender')

#px(g(y))Tracez la fonction convertie simplement comme
x = gaussianDist(sig,mu,g(y))
plt.plot(x/(x.sum()*0.01) ,y,'lime')

####
#Moyenne mu et g^-1(mu)Tracez la relation avec
plt.plot([mu, mu], [0, invg(mu)], 'k--')
plt.plot([0, mu], [invg(mu), invg(mu)], 'k--')

Résultat d'exécution