[PYTHON] Régression de processus gaussien Implémentation Numpy et GPy

introduction

Vous pouvez voir un graphique qui vous fait vous demander si vous pouvez tracer une telle ligne de régression. Par exemple

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Times New Roman"      #Définir la police entière
plt.rcParams["xtick.direction"] = "in"               #Vers l'intérieur de la ligne d'échelle de l'axe x
plt.rcParams["ytick.direction"] = "in"               #Vers l'intérieur de la ligne d'échelle de l'axe y
plt.rcParams["xtick.minor.visible"] = True           #Ajout de l'échelle auxiliaire sur l'axe x
plt.rcParams["ytick.minor.visible"] = True           #Ajout de l'échelle auxiliaire de l'axe y
plt.rcParams["xtick.major.width"] = 1.5              #Largeur de ligne de la ligne d'échelle principale de l'axe x
plt.rcParams["ytick.major.width"] = 1.5              #Largeur de ligne de la ligne d'échelle principale de l'axe y
plt.rcParams["xtick.minor.width"] = 1.0              #Largeur de ligne de la ligne d'échelle auxiliaire de l'axe x
plt.rcParams["ytick.minor.width"] = 1.0              #Largeur de ligne de la ligne d'échelle auxiliaire de l'axe y
plt.rcParams["xtick.major.size"] = 10                #longueur de la ligne d'échelle principale sur l'axe des x
plt.rcParams["ytick.major.size"] = 10                #Longueur de la ligne d'échelle principale de l'axe y
plt.rcParams["xtick.minor.size"] = 5                 #longueur de la ligne d'échelle auxiliaire sur l'axe des x
plt.rcParams["ytick.minor.size"] = 5                 #longueur de la ligne d'échelle auxiliaire sur l'axe y
plt.rcParams["font.size"] = 14                       #Taille de police
plt.rcParams["axes.linewidth"] = 1.5                 #Épaisseur du boîtier

def func(x):
    return 1.5*x+2.0 + np.random.normal(0.0,0.1,len(x))

x = np.array([0.0,0.2,0.5,1.0,10.0])
y = func(x) #Point de mesure
a,b = np.polyfit(x,y,1) #Régression linéaire

fig,axes = plt.subplots()
axes.scatter(x,y)
axes.plot(x,a*x+b)

J'ai pu tracer une ligne droite par régression linéaire. Cependant, il n'y a pas de points de mesure entre 1 et 10, et cette partie semble intuitivement incertaine. Lorsque j'ai recherché une bonne analyse de régression qui pourrait exprimer une telle incertitude, j'ai trouvé quelque chose appelé régression de processus gaussien, donc j'aimerais l'implémenter pour étude et l'essayer.

référence

Même moi, qui ne connais pas les mathématiques, pourrais le lire. Cet article correspond aux chapitres 2 et 3. Processus Gauss et apprentissage automatique

Les hyperparamètres existaient dans la régression des processus gaussiens et la méthode MCMC a été utilisée pour les estimer. Échantillonnage par la méthode Markov Chain Monte Carlo (MCMC) avec animateur

Il existe une bibliothèque sans avoir à l'implémenter vous-même. [Gpy vs scikit-learn: retour de processus gaussien avec python](https://nykergoto.hatenablog.jp/entry/2017/05/29/python%E3%81%A7%E3%82%AC%E3%82% A6% E3% 82% B9% E9% 81% 8E% E7% A8% 8B% E5% 9B% 9E% E5% B8% B0_ ~ % E3% 83% A2% E3% 82% B8% E3% 83% A5% E3% 83% BC% E3% 83% AB% E3% 81% AE% E6% AF% 94% E8% BC% 83 ~)

Régression du processus gaussien

1. Processus gaussien

Pour tout nombre naturel $ N $, le vecteur de sortie correspondant à l'entrée $ x_1, x_2, ..., x_N $ ${\mathbf f}=(f(x_1),f(x_2),...,f(x_N))$ Est la moyenne $ {\ mathbf \ mu} = (\ mu (x_1), \ mu (x_2), ..., \ mu (x_N)) $, $ K_ {nn ^ {'}} = k (x_n, x_n) ^ {'}) En suivant la distribution gaussienne $ {\ mathbf \ mu}, {\ rm K}) $ avec la matrice $ \ rm K $ comme élément et la matrice de covariance $ f On dit qu'il suit le processus gaussien $f \sim {\rm GP}(\mu({\mathbf x}), k({\mathbf x},{\mathbf x^{'}}))$ Écrire.

Cependant, je ne suis pas sûr de la nature du processus gaussien, donc je vais essayer de calculer la matrice du noyau en utilisant la fonction de base radiale comme fonction de noyau et un échantillon de la distribution gaussienne qui suit la matrice du noyau.

def rbf(x,s=[1.0,1.0,1.0e-6]):
    """
    Radial Basis Function, RBF :Fonction de base radiale
    """
    s1,s2,s3 = s[0],s[1],s[2]
    
    X, X_ = np.meshgrid(x,x)
    K = s1**2 * np.exp(-((X-X_)/s2)**2) + np.identity(len(x)) * s3**2
    return K

N = 100
x = np.linspace(0.0,10.0,N) #contribution
K = rbf(x) #Matrice du noyau

fig,axes = plt.subplots()
axes.imshow(K)

Pour échantillonner à partir d'une distribution gaussienne multidimensionnelle avec une matrice distribuée co-distribuée de $ \ rm K $, générez aléatoirement $ \ mathbf x $ et convertissez-le en $ \ mathbf {y} = \ mathbf {Lx} $ .. $ \ mathbf L $ peut être obtenu par une décomposition en clé collée de $ {\ rm K} $.

L = np.linalg.cholesky(K) #Décomposition Colleskey de la matrice du noyau

fig,axes = plt.subplots()
axes.plot()
for i in range(10):
    _x = np.random.normal(0.0,1.0,N)
    axes.plot(x,np.dot(L,_x))

Les fonctions apparaîtront au hasard, tout comme les nombres apparaissent lorsque vous lancez les dés. En outre, il semble que la nature intuitive du processus gaussien soit que des valeurs similaires proviennent de valeurs similaires.

2. Régression du processus gaussien

Il existe des mesures $ N . Par souci de simplicité, la valeur moyenne de y est 0. ${\mathcal D} = {(x_1,y_1),(x_2,y_2),...,(x_N,y_N)}$$ A ce moment, entre $ \ mathbf x $ et $ y $ $y=f({\mathbf x})$ Cette fonction $ f $ est un processus gaussien avec une moyenne de 0 $f \sim {\rm GP}(0,k(x,x^{'}))$ Supposons qu'il soit généré à partir de. Si vous définissez $ {\ mathbf y} = (y_1, y_2, ..., y_N) ^ {T} $, ce $ \ mathbf y $ suit une distribution gaussienne et la fonction noyau $ k $ pour chaque paire d'entrées. se servir de ${\rm K}\ =k(x,x^{'})$ En utilisant la matrice du noyau $ \ rm K $ donnée dans $y \sim {\mathcal N}(0,{\rm K})$ Est établi.

Voici comment trouver la valeur de $ y ^ {\ ast} $ à $ x ^ {\ ast} $ qui n'est pas incluse dans les données. $ \ Mathbf y $ plus $ y ^ {\ ast} $ $ {\ mathbf y} ^ {'} = (y_1, y_2, ..., y_N, y ^ {\ ast}) ^ { Soit T} $. Si vous faites $ \ rm K ^ {'} $ sur la matrice du noyau calculée à partir de $ (x_1, x_2, ..., x_N, x ^ {\ ast}) ^ {T} , tous suivent également la distribution gaussienne. Alors $y^{'} \sim {\mathcal N}(0,{\rm K}^{'})$ Ce sera. C'est, $\begin{pmatrix} {\mathbf y} \\ y^{\ast} \end{pmatrix} \sim {\mathcal N} \begin{pmatrix}0, {\begin{pmatrix}{\rm K} & {k_{\ast}} \\ {k_{\ast}^T} & {k_{\ast \ast}}\end{pmatrix}}\end{pmatrix}$$ Est établi. Puisqu'il s'agit d'une expression pour la distribution simultanée de $ y ^ {\ ast} $ et $ \ mathbf y $, la probabilité conditionnelle de $ y ^ {\ ast} $ étant donnée $ \ mathbf y $ est une distribution gaussienne. Il est calculé à partir de la probabilité conditionnelle entre les éléments de. La probabilité conditionnelle est: $p(y^{\ast}|x^{\ast},{\mathcal D}) = {\mathcal N}(k_{\ast}^T {\rm K}^{-1} {\mathbf y}, k_{\ast \ast}-k_{\ast}^T {\rm K}^{-1} k_{\ast})$

Écrivez le code selon la formule ci-dessus.

def rbf(x_train,s,x_pred=None):
    """
    Radial Basis Function, RBF :Fonction de base radiale
    """
    s1,s2,s3 = s[0],s[1],s[2]

    if x_pred is None:
        x = x_train
        X, X_ = np.meshgrid(x,x)
        K = s1**2 * np.exp(-((X-X_)/s2)**2) + np.identity(len(x)) * s3**2 #Matrice du noyau
        return K
    else:
        x = np.append(x_train,x_pred)
        X, X_ = np.meshgrid(x,x)
        K_all = s1**2 * np.exp(-((X-X_)/s2)**2) + np.identity(len(x)) * s3**2 #Matrice du noyau
        K = K_all[:len(x_train),:len(x_train)]
        K_s = K_all[:len(x_train),len(x_train):]
        K_ss = K_all[len(x_train):,len(x_train):]
        return K,K_s,K_ss
    
def pred(x_train,y_train,x_pred,s):
    K,K_s,K_ss = rbf(x_train,s,x_pred=x_pred)
    K_inv = np.linalg.inv(K) #Matrice inverse
    y_pred_mean = np.dot(np.dot(K_s.T,K_inv), y_train) #Valeur attendue de y
    y_pred_cov = K_ss - np.dot(np.dot(K_s.T,K_inv), K_s) #Matrice co-distribuée distribuée
    y_pred_std = np.sqrt(np.diag(y_pred_cov)) #écart-type
    return y_pred_mean,y_pred_std

Créez une fonction appropriée, échantillonnez-la et utilisez-la comme valeur mesurée.

def func(x):
    return np.sin(2.0*np.pi*0.2*x) + np.sin(2.0*np.pi*0.1*x)

pred_N = 100 #Nombre de points prédits
N = 10 #Nombre de points de mesure
x_train = np.random.uniform(0.0,10.0,N) #Point de mesure:x
y_train = func(x_train) + np.random.normal(0,0.1,1) #Point de mesure:y

x_true = x_pred = np.linspace(-2.0,12.0,pred_N)
fig,axes = plt.subplots()
axes.plot(x_true, func(x_true), label="True")
axes.scatter(x_train,y_train, label="Measured")
axes.legend(loc="best")

Faisons une régression de processus gaussien.

s = [1.0,1.0,0.1]
y_pred_mean,y_pred_std = pred(x_train,y_train,x_pred,s)

fig,axes = plt.subplots(figsize=(8,6))
axes.set_title("s1=1.0, s2=1.0, s3=0.1")
axes.plot(x_true, func(x_true), label="True")
axes.scatter(x_train,y_train, label="Measured")
axes.plot(x_pred, y_pred_mean, label="Predict")
plt.fill_between(x_pred,y_pred_mean+y_pred_std,y_pred_mean-y_pred_std,facecolor="b",alpha=0.3)
axes.legend(loc="best")
axes.set_xlim(-2.0,12.0)

J'ai pu faire un retour comme ça. La distribution est répandue là où le nombre de points de mesure est faible. Maintenant, les hyperparamètres de la fonction noyau sont donnés à la main. Voyons ce qui se passe si nous changeons cela.

Le résultat a beaucoup changé. Je veux aussi estimer les hyper paramètres.

Mettez les hyper paramètres ensemble comme $ \ theta $. La matrice du noyau dépend de $ \ theta $, où la probabilité de formation des données est

Si vous prenez le journal $\ln{p({\mathbf y}|{\mathbf x},\theta)} \propto -\ln{|{\rm K}|}-{\mathbf y}^T {\rm K}_{\theta}^{-1} {\mathbf y} + const.$ Et trouve $ \ theta $ qui maximise l'équation ci-dessus.

La méthode MCMC est utilisée car il est facile de se coincer dans la solution locale dans la méthode du gradient.

import emcee
def objective(s): #Fonction objective
    K = rbf(x_train,s)
    return -np.linalg.slogdet(K)[1] - y_train.T.dot(np.linalg.inv(K)).dot(y_train) #preuve

ndim = 3
nwalker = 6
s0 = np.random.uniform(0.0,5.0,[nwalker,ndim]) #Position initiale
sampler = emcee.EnsembleSampler(nwalker,ndim,objective) #Faire un échantillonneur
sampler.run_mcmc(s0,5000) #Commencer l'échantillonnage

s_dist = sampler.flatchain #Obtenez le résultat de l'échantillonnage
s = s_dist[sampler.flatlnprobability.argmax()]

y_pred_mean,y_pred_std = pred(x_train,y_train,x_pred,s)

fig,axes = plt.subplots(figsize=(8,6))
axes.plot(x_true, func(x_true), label="True")
axes.scatter(x_train,y_train, label="Measured")
axes.plot(x_pred, y_pred_mean, label="Predict")
plt.fill_between(x_pred,y_pred_mean+y_pred_std,y_pred_mean-y_pred_std,facecolor="b",alpha=0.3)
axes.legend(loc="best")
axes.set_xlim(-2.0,12.0)

Ça fait du bien.

3. Régression de processus gaussien à l'aide de la bibliothèque GPy

Il existe une bibliothèque pratique de régression de processus gaussien. Il est recommandé car il possède une grande variété de noyaux et est facile à visualiser.

import GPy
import GPy.kern as gp_kern

kern = gp_kern.RBF(input_dim=1)
gpy_model = GPy.models.GPRegression(X=x_train.reshape(-1, 1), Y=y_train.reshape(-1, 1), kernel=kern, normalizer=None)

fig,axes = plt.subplots(figsize=(8,6))
axes.plot(x_true, func(x_true), c="k", label="True")
gpy_model.optimize()
gpy_model.plot(ax=axes)
axes.legend(loc="best")
axes.set_xlim(-2.0,12.0)

Résumé

La régression du processus gaussien a permis d'estimer l'incertitude comprise. Dans cet article, je n'ai pu considérer que la partie théorie de base, mais la régression du processus gaussien a un problème avec la quantité de calcul, et il y a diverses discussions sur la façon de le résoudre. De plus, il existe de nombreux types de fonctions du noyau, pas seulement un type, et il est nécessaire de les sélectionner et de les combiner de manière appropriée en fonction de la cible d'analyse.