[PYTHON] Normalisation de la théorie et de la mise en œuvre des flux

introduction

J'ai implémenté la normalisation du flux, qui est l'une des méthodes d'inférence de variables, je vais donc la laisser ici sous forme de mémorandum. Nous commencerons par la régression linéaire bayésienne afin de pouvoir couvrir tout le flux, puis nous expliquerons étape par étape pourquoi l'inférence de variable est nécessaire et comment la normalisation du flux effectue l'inférence de variable. Bien qu'il s'agisse d'un flux de normalisation, il peut également être utilisé comme modèle de génération en plus de l'inférence de variables. GAN et VAE sont connus comme des modèles de génération, mais la normalisation du flux appartient à un modèle de génération différent basé sur le flux. L'article sur la normalisation du flux est ici. Le code implémenté est publié sur GitHub.

Régression linéaire bayésienne

Commençons par l'histoire de la régression linéaire bayésienne. Pour récapituler, dans la régression linéaire normale, la valeur prédite $ \ boldsymbol y $ est sortie pour les données observées $ \ boldsymbol x $ comme indiqué ci-dessous. $ \ Boldsymbol \ theta $ est un paramètre et correspond à la pente d'une droite dans le cas d'une régression simple. $\boldsymbol y=\boldsymbol \theta^T\boldsymbol x \tag{1}$ Cependant, comme la valeur réelle contient du bruit, pour modéliser une telle ambiguïté, nous supposons une distribution gaussienne avec une moyenne de $ \ boldsymbol y $ et une covariance de $ I $ comme suit. (Par exemple, lors de la mesure de la température, la valeur obtenue par l'erreur d'équipement du thermomètre comprend le bruit, c'est donc une image qui exprime le bruit par la dispersion de la distribution gaussienne.) $p(\boldsymbol y \mid \boldsymbol x，\boldsymbol \theta)=\mathcal{N}(\boldsymbol y\mid \boldsymbol \theta^T\boldsymbol x， I) \tag{2}$ Ici, afin de faire une prédiction de régression pour les données non observées en utilisant l'équation ci-dessus, il est nécessaire de déterminer le paramètre optimal $ \ boldsymbol \ theta $. À partir de là, nous adoptons une approche bayésienne, mais Bayes traite le paramètre $ \ boldsymbol \ theta $ comme une distribution de probabilité. Autrement dit, la distribution de probabilité $ p (\ boldsymbol \ theta) du paramètre optimal $ \ boldsymbol \ theta $ à partir des données $ \ boldsymbol X ^ {(train)} $ et $ \ boldsymbol y ^ {(train)} $ déjà à portée de main \ mid \ boldsymbol X ^ {(train)}, \ boldsymbol y ^ {(train)}) Détermine $. (L'approche non-Bayes trouve le paramètre optimal $ \ boldsymbol \ theta $ comme une valeur déterministe plutôt qu'une distribution de probabilité. Dans ce cas, une technique appelée estimation la plus probable est utilisée, ce qui donne une méthode des moindres carrés. ) $ P (\ boldsymbol \ theta \ mid \ boldsymbol X ^ {(train)}, \ boldsymbol y ^ {(train)}) $ est calculé en utilisant le théorème bayésien suivant. $ p(\boldsymbol \theta \mid \boldsymbol X^{(train)}，\boldsymbol y^{(train)})= \cfrac{p(\boldsymbol y^{(train)} \mid \boldsymbol X^{(train)}， \boldsymbol \theta)p(\boldsymbol \theta)}{p(\boldsymbol y^{(train)} \mid \boldsymbol X^{(train)})} \tag{3} $ $ p (\ boldsymbol \ theta \ mid \ boldsymbol X ^ {(train)}, \ boldsymbol y ^ {(train)}) $ est appelée la distribution postérieure, et une fois la distribution postérieure obtenue, la formule suivante est utilisée. Vous pouvez l'utiliser pour calculer la prédiction de régression $ \ boldsymbol y ^ {(nouveau)} $ pour les données inconnues $ \ boldsymbol X ^ {(nouveau)} . (Bien que cela semble compliqué, je n'utiliserai pas ces formules directement cette fois, alors veuillez la lire à la légère.) $ p(\boldsymbol y^{(new)}\mid \boldsymbol X^{(new)}， \boldsymbol y^{(train)}， \boldsymbol X^{(train)})= \int p(\boldsymbol y^{(new)}\mid \boldsymbol X^{(new)}，\boldsymbol \theta) p(\boldsymbol \theta \mid \boldsymbol y^{(train)}， \boldsymbol X^{(train)}) d\boldsymbol \theta $$ Revenant à l'histoire, regardons de plus près le théorème de Bayes (3) afin d'obtenir la distribution postérieure. Premièrement, la molécule $ p (\ boldsymbol y ^ {(train)} \ mid \ boldsymbol X ^ {(train)}, \ boldsymbol \ theta) $ est appelée vraisemblance et a déjà été modélisée par l'équation (2). Comme c'est déjà fait, il peut être calculé. De même, $ p (\ boldsymbol \ theta) $ apparaissant dans la molécule de l'équation (3) est appelé une distribution a priori, qui est modélisée par le concepteur sur la base de la connaissance préalable de la distribution de $ \ boldsymbol \ theta $. Je vais. (Il existe des règles pour choisir, mais ce n'est pas une erreur de choisir une distribution de probabilité telle qu'une distribution gaussienne ou une distribution uniforme.)

Enfin, $ p (\ boldsymbol y ^ {(train)} \ mid \ boldsymbol X ^ {(train)}) , qui est appelé preuve ou probabilité marginale, est calculé à partir de l'équation suivante. $ p(\boldsymbol y^{(train)} \mid \boldsymbol X^{(train)})= \int p(\boldsymbol y^{(train)} \mid \boldsymbol X^{(train)}， \boldsymbol \theta)p(\boldsymbol \theta)d\boldsymbol \theta \tag{4} $$ Il est possible de calculer cette équation de manière analytique avec un modèle simple tel que la régression linéaire, mais la complexité du modèle rend difficile le calcul de l'intégrale. La méthode utilisée est l'inférence différentielle. (Il existe d'autres méthodes d'approximation telles que la méthode de Monte Carlo.)

Raisonnement variant

En évitant l'analyse analytique exacte de la vraisemblance marginale, essayez d'approximer directement la distribution postérieure $ p (\ boldsymbol z \ mid \ boldsymbol x) $ que vous voulez trouver par la distribution de probabilité $ q (\ boldsymbol z) $. L'idée est un raisonnement variable. (Afin de faire correspondre la notation avec le papier, le paramètre $ \ boldsymbol \ theta $ qui est apparu jusqu'à présent est décrit comme $ \ boldsymbol z $.) Comme exemple d'inférence de variable, il existe une méthode de décomposition de la distribution de probabilité en supposant l'indépendance de la distribution pour la distribution postérieure approximative $ q (\ boldsymbol z) . Cette méthode est appelée approximation du champ moyen et prend la forme d'une distribution postérieure approximative comme le montre l'équation suivante. $q(\boldsymbol z)=\prod_{i} q(\boldsymbol z_i) \tag{5}$$ Ensuite, une optimisation est effectuée pour rendre cette distribution de probabilité $ q (\ boldsymbol z) $ proche de la vraie distribution postérieure $ p (\ boldsymbol z \ mid \ boldsymbol x) $. L'optimisation de l'inférence différentielle implique une transformation d'expression légèrement difficile, donc pour conclure d'abord, maximisez la fonction de perte de ELBO (Evidence Lower Bound) $ \ mathcal {L} (\ boldsymbol {x}) $ dans l'équation suivante. Il optimise en faisant. $ \mathcal{L}(\boldsymbol{x})=\int q(\boldsymbol z)\ln \frac{p(\boldsymbol x， \boldsymbol z)}{q(\boldsymbol z)}d\boldsymbol z \tag{6} $ De là, ELBO sera dérivé, mais vous pouvez l'ignorer. Premièrement, le but initial est de trouver une distribution postérieure approximative $ q (\ boldsymbol z) $ qui maximise la preuve $ p (\ boldsymbol x) $, qui est la probabilité de génération de données. Maximiser $ p (\ boldsymbol x) $ équivaut à maximiser $ \ ln p (\ boldsymbol x) $, qui est son logarithme, alors considérez ce dernier pour plus de simplicité. La formule suivante est obtenue en procédant au calcul pour $ \ ln p (\ boldsymbol x) $.

\begin{eqnarray}
\ln p(\boldsymbol x)=\ln \int p(\boldsymbol x， \boldsymbol z)d\boldsymbol z\\
=\ln \int q(\boldsymbol z)\frac{p(\boldsymbol x， \boldsymbol z)}{q(\boldsymbol z)}d\boldsymbol z\\
\geq \int q(\boldsymbol z)\ln \frac{p(\boldsymbol x， \boldsymbol z)}{q(\boldsymbol z)}d\boldsymbol z \\
= \mathcal{L}(\boldsymbol{x})
\end{eqnarray}

Voici la transformation de la formule à l'étape (3), mais elle a été calculée en utilisant l'inégalité de Jensen. Comme vous pouvez le voir dans cette formule, la borne inférieure de $ \ ln p (\ boldsymbol x) $ est ELBO $ \ mathcal {L} (\ boldsymbol {x}) $. Je n'entrerai pas dans les détails, mais maximiser ELBO équivaut à minimiser la divergence KL des distributions postérieures et postérieures approximatives cibles, donc cette tâche peut être réduite à la maximisation ELBO.

Revenant à l'histoire, supposons que la distribution postérieure approximative $ q (\ boldsymbol z) $ puisse être obtenue en résolvant le problème d'optimisation. Cependant, en premier lieu, il y a un problème avec l'hypothèse de l'équation (5). En supposant l'indépendance comme dans l'équation (5), la corrélation existant entre les variables ne peut pas être considérée, et la distribution de probabilité devient peu expressive, et la vraie distribution postérieure ne peut pas être exprimée. .. La normalisation du flux est l'une des méthodes d'inférence variable qui résout ce problème de manque d'expressivité.

Normalizing Flow Normaliser Flow est une distribution complexe $ q_k (\ boldsymbol z_k) en superposant une transformation non linéaire $ f $ sur une variable de probabilité $ \ boldsymbol z $ qui suit une distribution de probabilité simple (distribution gaussienne, etc.) $ q (\ boldsymbol z) $. ) L'idée est d'obtenir $. <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/302831/4cfcd80e-18ee-f050-df70-5f59a8e3e8fd.png ", width="100%"> Source: L.Weng. "Modèle génératif profond basé sur le flux

Regardons de plus près. Les équations (7) et (8) sont obtenues en appliquant la transformation non linéaire $ f $ à la variable de probabilité $ \ boldsymbol z $ qui suit la distribution de probabilité $ q (\ boldsymbol z) . Le différentiel qui apparaît sur le côté droit ici est jacobien. det signifie une expression matricielle. Cependant, à la fin, vous pouvez calculer le Jacobien et le Det manuellement, vous n'avez donc pas à réfléchir sérieusement. $\boldsymbol z'=f(\boldsymbol z)\tag{7}$ $q'(\boldsymbol z')=q(\boldsymbol z)\left| \det \frac{\partial f}{\partial \boldsymbol z'}\right|^{-1} \tag{8}$ En répétant cette transformation non linéaire plusieurs fois, la distribution de probabilité devient plus complexe et expressive. $ \boldsymbol z_k=f_k\circ \cdots \circ f_2 \circ f_1(\boldsymbol z_0) \tag{9} $ $ q_k(\boldsymbol z_k)=q_0(\boldsymbol z_0)\prod_{k=1}^{K} \left| \det \frac{\partial f_k}{\partial \boldsymbol z_{k-1}}\right|^{-1} \tag{10} $$ En ce qui concerne la sélection de la fonction non linéaire $ f_k $, c'est OK s'il s'agit d'une fonction entièrement à un seul coup, c'est-à-dire une fonction non linéaire qui peut calculer la fonction inverse. Dans l'article, deux types de fonction non linéaire $ f_k $ sont introduits: l'écoulement planaire et l'écoulement radial. Ici, nous n'expliquerons que l'écoulement planaire.

l'écoulement planaire est la fonction suivante. $ \ boldsymbol u, \ boldsymbol w et b $ sont des paramètres et leurs valeurs sont mises à jour par apprentissage. $ h $ est une fonction d'activation et utilise tanh. $ f(\boldsymbol z)=\boldsymbol z + \boldsymbol u h(\boldsymbol w^T \boldsymbol z + b) \tag{11} $ En utilisant la fonction de l'équation ci-dessus, l'équation de la matrice jacobienne peut être calculée comme suit. $ \left| \det \frac{\partial f(\boldsymbol z)}{\partial \boldsymbol z}\right| =\left| 1+\boldsymbol u^T \psi(\boldsymbol z)\right| \tag{12} $ Cependant, $ \ psi (\ boldsymbol z) $ est calculé à partir de l'équation suivante. $ \psi(\boldsymbol z)=h'(\boldsymbol w^T \boldsymbol z+b)\boldsymbol w \tag{13} $ Cela permet de calculer l'équation (10), et la distribution de probabilité après conversion peut être obtenue. $ q_k(\boldsymbol z_k)=q_0(\boldsymbol z_0)\prod_{k=1}^{K} \left| 1+\boldsymbol u_k^T \psi_k(\boldsymbol z_k)\right|^{-1} \tag{13} $ L'équation (13) est implémentée en écoulement planaire, mais vous pouvez voir que le coût de calcul est très faible (il suffit de calculer le produit intérieur des vecteurs et de prendre l'inverse de la valeur absolue).

Ensuite, nous optimisons la distribution de probabilité obtenue par l'équation (13) pour la rapprocher de la vraie distribution postérieure. Comme mentionné précédemment, dans l'inférence différentielle, l'optimisation est effectuée en maximisant ELBO. Nous calculerons manuellement la formule ELBO pour normaliser le débit.

\mathcal{L}(\boldsymbol{x})=
\int q(\boldsymbol z)\ln \frac{p(\boldsymbol x， \boldsymbol z)}{q(\boldsymbol z)}d\boldsymbol z\\
=\mathbb{E}_{q(\boldsymbol z)}[\ln p(\boldsymbol x， \boldsymbol z) - \ln q(\boldsymbol z)]\\
=\mathbb{E}_{q_k(\boldsymbol z_k)}[\ln p(\boldsymbol x， \boldsymbol z_k) - \ln q(\boldsymbol z_k)]\\
\sim \sum_{l=1}^{L}[\ln p(\boldsymbol{x}， \boldsymbol{z_k^{(l)}})-\ln q(\boldsymbol{z_k^{(l)}})]
\tag{14}

La somme qui apparaît à la fin de l'expression est la somme des L exemples de variables probabilistes contenues dans le mini-lot (indice $ l). $ Signifie le $ l $ ème échantillon). Lors de la mise en œuvre, la minimisation est plus pratique que la maximisation, de sorte que l'équation suivante, qui est ELBO multipliée par un moins, est minimisée. $ -\mathcal{L}(\boldsymbol{x})=\sum_{l=1}^{L}[\ln q(\boldsymbol{z_k^{(l)}}) - \ln p(\boldsymbol{x}， \boldsymbol{z_k^{(l)}})]　\tag{15} $ C'est la fin de l'explication de la théorie. Le reste est la mise en œuvre.

Environnement d'exécution

Python 3.7.5 TensorFlow 1.15.0

la mise en oeuvre

Implémenter par TensorFlow. Comme pour l'article, le problème est la distribution gaussienne comme distribution initiale, et la distribution cible est reproduite en superposant des transformations non linéaires par écoulement planaire. (La gauche est la distribution gaussienne, qui est la distribution initiale, et le centre et la droite sont les distributions cibles).

Tout d'abord, nous allons l'implémenter à partir du flux planaire. Puisque $ \ ln q_k (\ boldsymbol z_k) $ est requis lors du calcul de la fonction de perte, il est calculé ici.

normalizing_flow.py

class PlanarFlow:
    def __init__(self, dim):
        self.dim = dim
        self.h = lambda x: tf.tanh(x)
        self.h_prime = lambda x: 1 - tf.tanh(x)**2
        self.w = tf.Variable(tf.random.truncated_normal(shape=(1, self.dim)))
        self.b = tf.Variable(tf.zeros(shape=(1)))
        self.u = tf.Variable(tf.random.truncated_normal(shape=(1, self.dim)))
        

    def __call__(self, z, log_q):
        z = z + self.u*self.h(tf.expand_dims(tf.reduce_sum(z*self.w, -1), -1) + self.b)
        psi = self.h_prime(tf.expand_dims(tf.reduce_sum(z*self.w, -1), -1) + self.b)*self.w
        det_jacob = tf.abs(1 + tf.reduce_sum(psi*self.u, -1))
        log_q = log_q -  tf.log(1e-7 + det_jacob)
        return z, log_q

Ensuite, un flux de normalisation est construit en empilant K flux planaires.

normalizing_flow.py

class NormalizingFlow:
    def __init__(self, K, dim):
        self.K = K
        self.dim = dim
        self.planar_flow = [PlanarFlow(self.dim) for i in range(self.K)]
        

    def __call__(self, z_0,log_q_0):
        z, log_q = self.planar_flow[0](z_0,log_q_0)
        for pf in self.planar_flow[1:]:
            z, log_q = pf(z, log_q)
        return z, log_q

Ensuite, la fonction de perte (15) est calculée comme suit. Ici, target_density calcule la distribution de probabilité cible. La méthode d'optimisation est également décidée ici. Utilisez Adam pour effectuer l'optimisation afin de minimiser la fonction de perte. Définissez également une fonction pour obtenir l'espace réservé.

normalizing_flow.py

def calc_loss(z_k, log_q_k, target_density):
    log_p = tf.log(target_density.calc_prob_tf(z_k)+1e-7)
    loss = tf.reduce_mean(log_q_k - log_p, -1)
    return loss


def get_train(loss):
    return tf.train.AdamOptimizer().minimize(loss)

def get_placeholder():
    z_0 = tf.placeholder(tf.float32, shape=[None, 2])
    log_q_0 = tf.placeholder(tf.float32, shape=[None])
    return z_0,log_q_0

Construisez un graphe de calcul en utilisant les classes et fonctions ci-dessus.

main.py

normalizing_flow = NormalizingFlow(K=16, dim=2)

z_0,log_q_0 = get_placeholder()
z_k, log_q_k = normalizing_flow(z_0,log_q_0)
loss = calc_loss(z_k, log_q_k, target_density)
train = get_train(loss)

L'apprentissage est effectué en utilisant la méthode de descente de gradient stochastique. Le mini-lot contient 1000 échantillons et le nombre d'apprentissage est de 100 000.

main.py

with tf.Session() as sess:
    invisible_axis = True
    sess.run(tf.global_variables_initializer())

    for iteration in range(100000+1):
        z_0_batch = normal_distribution.sample(1000)
        log_q_0_batch = np.log(normal_distribution.calc_prob(z_0_batch))
        _, loss_value = sess.run([train, loss], {z_0:z_0_batch, log_q_0:log_q_0_batch})

Le résultat est là. Vous pouvez voir que la distribution complexe est échantillonnée correctement, à partir de la distribution gaussienne. <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/302831/4ab8c588-1d49-e110-c3c9-deb72f6e8eac.png "， width=33%> <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/302831/4d21dd7d-6aa1-fd70-57c0-e689f7525111.png "， width=33%><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/302831/4dc89cef-d8ac-a2b6-4897-4224d1de1fae.png "， width=33%>

Toutes les implémentations

Enfin, voici un résumé de toutes les implémentations. Il existe trois types: main.py, normalizing_flow.py et distribution.py. main.py est un fichier pour construire un graphe de calcul et exécuter un entraînement, normalizing_flow.py est un fichier pour définir un modèle et une fonction de perte pour construire un graphe de calcul, et enfin distribution.py est un échantillon de distribution de probabilité. C'est un fichier qui a une fonction pour calculer la probabilité. Le fichier est également publié sur GitHub.

main.py

import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt
from distribution import *
from normalizing_flow import *

normal_distribution = NormalDistribution2D()
target_density = TargetDistribution1()

normalizing_flow = NormalizingFlow(K=16, dim=2)

z_0,log_q_0 = get_placeholder()
z_k, log_q_k = normalizing_flow(z_0,log_q_0)
loss = calc_loss(z_k, log_q_k, target_density)
train = get_train(loss)

with tf.Session() as sess:
    invisible_axis = True
    sess.run(tf.global_variables_initializer())

    for iteration in range(100000+1):
        z_0_batch = normal_distribution.sample(1000)
        log_q_0_batch = np.log(normal_distribution.calc_prob(z_0_batch))
        _, loss_value = sess.run([train, loss], {z_0:z_0_batch, log_q_0:log_q_0_batch})
        
        if iteration % 100 == 0:
            print('Iteration : {}   Loss : {}'.format(iteration, loss_value))
            
        if iteration % 10000 == 0:
            z_k_value = sess.run(z_k, {z_0:z_0_batch, log_q_0:log_q_0_batch})
            plt.figure(figsize=(6, 6))
            plt.scatter(z_k_value[:, 0], z_k_value[:, 1], alpha=0.7)
            if invisible_axis:
                plt.tick_params(bottom=False,left=False,right=False,top=False)
                plt.tick_params(labelbottom=False,labelleft=False,labelright=False,labeltop=False)
            plt.show()

normalizing_flow.py

import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt


def get_placeholder():
    z_0 = tf.placeholder(tf.float32, shape=[None, 2])
    log_q_0 = tf.placeholder(tf.float32, shape=[None])
    return z_0,log_q_0


def calc_loss(z_k, log_q_k, target_density):
    log_p = tf.log(target_density.calc_prob_tf(z_k)+1e-7)
    loss = tf.reduce_mean(log_q_k - log_p, -1)
    return loss


def get_train(loss):
    return tf.train.AdamOptimizer().minimize(loss)


class PlanarFlow:
    def __init__(self, dim):
        self.dim = dim
        self.h = lambda x: tf.tanh(x)
        self.h_prime = lambda x: 1 - tf.tanh(x)**2
        self.w = tf.Variable(tf.random.truncated_normal(shape=(1, self.dim)))
        self.b = tf.Variable(tf.zeros(shape=(1)))
        self.u = tf.Variable(tf.random.truncated_normal(shape=(1, self.dim)))
        

    def __call__(self, z, log_q):
        z = z + self.u*self.h(tf.expand_dims(tf.reduce_sum(z*self.w, -1), -1) + self.b)
        psi = self.h_prime(tf.expand_dims(tf.reduce_sum(z*self.w, -1), -1) + self.b)*self.w
        det_jacob = tf.abs(1 + tf.reduce_sum(psi*self.u, -1))
        log_q = log_q -  tf.log(1e-7 + det_jacob)
        return z, log_q


class NormalizingFlow:
    def __init__(self, K, dim):
        self.K = K
        self.dim = dim
        self.planar_flow = [PlanarFlow(self.dim) for i in range(self.K)]
        
        
    def __call__(self, z_0,log_q_0):
        z, log_q = self.planar_flow[0](z_0,log_q_0)
        for pf in self.planar_flow[1:]:
            z, log_q = pf(z, log_q)
        return z, log_q

distribution.py

import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt

class Distribution:    
    def calc_prob(self, z):
        p = np.zeros(z.shape[0])
        return p
    

    def plot(self, size=5):
        side = np.linspace(-size, size, 1000)
        z1, z2 = np.meshgrid(side, side)
        shape = z1.shape
        z1 = z1.ravel()
        z2 = z2.ravel()
        z = np.c_[z1, z2]
        probability = self.calc_prob(z).reshape(shape)
        plt.figure(figsize=(6, 6))
        plt.imshow(probability)
        plt.tick_params(bottom=False,left=False,right=False,top=False)
        plt.tick_params(labelbottom=False,labelleft=False,labelright=False,labeltop=False)
        plt.show()


class NormalDistribution2D(Distribution):
    def sample(self, sample_num):
        z = np.random.randn(sample_num, 2)
        return z
    

    def sample_tf(self, sample_num):
        z = tf.random_normal([sample_num, 2])
        return z
    

    def calc_prob(self, z):
        p = np.exp(-(z[:, 0]**2+z[:, 1]**2)/2)/(2*np.pi) 
        return p
    

    def calc_prob_tf(self, z):
        p = tf.exp(-(z[:, 0]**2+z[:, 1]**2)/2)/(2*np.pi) 
        return p


class TargetDistribution1(Distribution):
    def calc_prob(self, z):
        z1, z2 = z[:, 0], z[:, 1]
        norm = np.sqrt(z1**2+z2**2)
        exp1 = np.exp(-0.5*((z1-2)/0.6)**2)
        exp2 = np.exp(-0.5*((z1+2)/0.6)**2)
        p = 0.5*((norm - 2)/0.4)**2 - np.log(exp1 + exp2)
        return np.exp(-p)
    

    def calc_prob_tf(self, z):
        z1, z2 = z[:, 0], z[:, 1]
        norm = tf.sqrt(z1**2+z2**2)
        exp1 = tf.exp(-0.5*((z1-2)/0.6)**2)
        exp2 = tf.exp(-0.5*((z1+2)/0.6)**2)
        p = 0.5*((norm - 2)/0.4)**2 - tf.log(exp1 + exp2)
        return tf.exp(-p)


class TargetDistribution2(Distribution):
    def calc_prob(self, z):
        z1, z2 = z[:, 0], z[:, 1]
        w1 = np.sin(0.5*np.pi*z1)
        p = 0.5*((z2 - w1)/0.4)**2
        return np.exp(-p)
    
    
    def calc_prob_tf(self, z):
        z1, z2 = z[:, 0], z[:, 1]
        w1 = tf.sin(0.5*np.pi*z1)
        p = 0.5*((z2 - w1)/0.4)**2
        return tf.exp(-p)

Les références

[1] Ian.Goodfellow, et.al., "Deep learning". [2] C.M. Bishop, "Reconnaissance de formes et apprentissage automatique". [3] C.M. Bishop, "Reconnaissance de formes et apprentissage automatique". [4] Atsushi Suyama, "Introduction à l'apprentissage automatique par l'inférence bayésienne".