[PYTHON] Implémenter le processus gaussien dans Pyro

Le processus gaussien était à l'origine l'un des processus stochastiques en temps continu, mais il a été appliqué à l'apprentissage automatique en tant que modèle probabiliste. Le fait est que ** les processus gaussiens sont des distributions gaussiennes multivariées de dimension infinie **. S'il s'agit d'une fonction, chaque sortie pour différentes entrées suit une distribution gaussienne. Le noyau définit la similitude des entrées. Si les entrées sont similaires, les sorties seront similaires, vous pouvez donc prédire une courbe lisse. D'un autre point de vue, vous pouvez le considérer comme un modèle de régression linéaire dans lequel les fonctions de base sont disposées à l'infini sur la grille. Compte tenu des données observées, la matrice du noyau est calculée pour déterminer la moyenne et la variance des emplacements non observés. (La distribution gaussienne ne change pas.) Dans le calcul, vous pouvez obtenir la sortie sans calculer le vecteur de caractéristiques par l'astuce du noyau.

Pour une explication théorique du processus gaussien, reportez-vous à «Processus gaussien et apprentissage automatique (série professionnelle d'apprentissage automatique)». Fait intéressant, il s'avère que le processus gaussien équivaut à un réseau de neurones avec un nombre infini d'unités (https://oumpy.github.io/blog/2020/04/neural_tangents.html).

Cette fois, je vais implémenter le processus gaussien avec Pyro présenté dans Article précédent.

Processus gaussien retour avec Pyro

Nous ferons Exemples du tutoriel officiel de Pyro. Le Document officiel est également utile.

import matplotlib.pyplot as plt
import torch
import pyro
import pyro.contrib.gp as gp
import pyro.distributions as dist
pyro.set_rng_seed(100)

Prenez 20 points de y = 0,5 * sin (3x) et utilisez-les comme données d'observation.

N = 20
X = dist.Uniform(0.0, 5.0).sample(sample_shape=(N,))
y = 0.5 * torch.sin(3*X) + dist.Normal(0.0, 0.2).sample(sample_shape=(N,))
plt.plot(X.numpy(), y.numpy(), 'kx')

À partir de ces 20 données d'observation, nous prédisons la fonction d'origine (dans ce cas, y = 0,5 * sin (3x)). Choisissez un noyau et effectuez une régression de processus gaussien.

#Définir les hyper paramètres
variance = torch.tensor(0.1)
lengthscale = torch.tensor(0.1)
noise = torch.tensor(0.01)
#Revenir
kernel = gp.kernels.RBF(input_dim=1, variance=variance, lengthscale=lengthscale)
gpr = gp.models.GPRegression(X, y, kernel, noise=noise)

Maintenant vous avez un retour. Affichez le résultat de la prédiction.

Xtest = torch.linspace(-0.5, 5.5, 500)
with torch.no_grad():
    mean, cov = gpr(Xtest, full_cov=True, noiseless=False)
sd = cov.diag().sqrt()
plt.plot(Xtest.numpy(), mean.numpy(), 'r', lw=2)
plt.fill_between(Xtest.numpy(), (mean - 2.0 * sd).numpy(), (mean + 2.0 * sd).numpy(), color='C0', alpha=0.3)
plt.plot(X.numpy(), y.numpy(), 'kx')

Dans la régression de processus gaussien, la fonction de prédiction est exprimée sous la forme d'un nuage de fonctions (un ensemble de courbes de prédiction) comme le montre la figure ci-dessus. Ce point est similaire au résultat de la régression linéaire bayésienne. Cependant, ce résultat est le résultat lorsque des hyper paramètres (variance, échelle de longueur, bruit) tels que le noyau sont déterminés de manière appropriée. Optimisez avec la méthode de descente de gradient de l'inférence de variantes pour ajuster les hyperparamètres de noyau appropriés. (MCMC est également acceptable)

optimizer = torch.optim.Adam(gpr.parameters(), lr=0.005)
loss_fn = pyro.infer.Trace_ELBO().differentiable_loss
losses = []
num_steps = 2500
for i in range(num_steps):
    optimizer.zero_grad()
    loss = loss_fn(gpr.model, gpr.guide)
    loss.backward()
    optimizer.step()
    losses.append(loss.item())
#Tracer la courbe de perte
plt.plot(losses)
plt.xlabel("step")
plt.ylabel("loss")
#Afficher les résultats d'optimisation
print('variance = {}'.format(gpr.kernel.variance))
print('lengthscale = {}'.format(gpr.kernel.lengthscale))
print('noise = {}'.format(gpr.noise))

variance = 0.15705525875091553 lengthscale = 0.4686208963394165 noise = 0.017524730414152145

Affichez le résultat de la prédiction.

Xtest = torch.linspace(-0.5, 5.5, 500)
with torch.no_grad():
    mean, cov = gpr(Xtest, full_cov=True, noiseless=False)
sd = cov.diag().sqrt()
plt.plot(Xtest.numpy(), mean.numpy(), 'r', lw=2)
plt.fill_between(Xtest.numpy(), (mean - 2.0 * sd).numpy(), (mean + 2.0 * sd).numpy(), color='C0', alpha=0.3)
plt.plot(X.numpy(), y.numpy(), 'kx')

La forme du nuage est plus lisse que dans la figure précédente, et elle est plus proche de la fonction d'origine (y = 0,5 * sin (3x)). Je pense que le résultat de la prédiction a été amélioré en ajustant les hyper paramètres.

Autres applications

Optimisation bayésienne

https://pyro.ai/examples/bo.html Vous pouvez trouver la valeur minimale lors de la recherche en émettant une fonction prédictive par régression (méthode de solution approximative).

Modèle de variable latente de processus gaussien (GPLVM)

https://pyro.ai/examples/gplvm.html GPLVM est un apprentissage non supervisé et peut réduire les dimensions en recherchant l'entrée à partir de la sortie.

À la fin

Dans l'analyse normale des données, vous devez décider du modèle probabiliste qui convient à votre situation. Par exemple, la modélisation statistique bayésienne définit un modèle pour déterminer la distribution antérieure. À cet égard, le processus gaussien n'a besoin que de spécifier le ** noyau (similarité entre les échantillons) **, il semble donc être très polyvalent. Et même ce noyau semble être automatiquement déterminé par ARD (détermination automatique de la pertinence).

De plus, le processus gaussien peut exprimer l'incertitude, et je pense qu'il est possible d'apprendre de manière flexible même avec un petit nombre d'échantillons. Cependant, c'est un modèle compliqué et difficile à interpréter.