Aperçu

Bonjour, c'est kwashi. Ces dernières années, avec la popularité de l'apprentissage automatique, il est devenu appliqué au domaine industriel, comme l'utilisation de la prédiction à l'aide d'un modèle de régression pour le contrôle. Cependant, l'incertitude prédictive est également importante dans les contextes critiques des produits. Par exemple, supposons que vous souhaitiez renvoyer les paramètres d'un appareil de production à partir de l'environnement d'usine (température, etc.). Cependant, même si un modèle qui n'est pas utilisé pour l'entraînement se produit soudainement, certains paramètres seront renvoyés. Par conséquent, l'incertitude de l'inférence est également importante quant à la fiabilité de la régression.

La statistique bayésienne est un domaine qui montre cette incertitude sous forme de distribution de probabilité. Les statistiques bayésiennes sont utilisées dans divers domaines tels que la conduite automatisée et le traitement acoustique. Cet article fournit un exemple de régression utilisant les statistiques bayésiennes. Retirez la fonction sin, qui est un modèle simple comme cible. Tout d'abord, un exemple de régression polypoly est montré. Cependant, en polypoly, la fonction sin ne peut pas être exprimée sans définir suffisamment de paramètres. Par conséquent, nous développerons une méthode pour estimer la fonction elle-même en utilisant le processus gaussien.

En outre, cet article se concentre sur les bibliothèques PyMC et Scipy de python.

Cet article ne couvre pas suffisamment les statistiques bayésiennes. J'ai écrit cet article en référence aux livres suivants, veuillez donc le lire. [Suite / Reconnaissance de formes facile à comprendre - Introduction à l'apprentissage sans enseignant](https://www.amazon.co.jp/%E7%B6%9A%E3%83%BB%E3%82%8F%E3%81%8B% E3% 82% 8A% E3% 82% 84% E3% 81% 99% E3% 81% 84% E3% 83% 91% E3% 82% BF% E3% 83% BC% E3% 83% B3% E8% AA% 8D% E8% AD% 98% E2% 80% 95% E6% 95% 99% E5% B8% AB% E3% 81% AA% E3% 81% 97% E5% AD% A6% E7% BF% 92% E5% 85% A5% E9% 96% 80% E2% 80% 95-% E7% 9F% B3% E4% BA% 95-% E5% 81% A5% E4% B8% 80% E9% 83% 8E / dp / 427421530X) Bayesian Analysis with Python

Régression polymérique de la fonction Sin

Dans ce chapitre, les données sont générées selon une distribution normale avec un écart-type moyen sin (x) de 0,2 à un angle x extrait aléatoirement d'une distribution uniforme. Il s'intègre également dans un polypole de cinquième ordre en utilisant la méthode des moindres carrés. Le polynôme du cinquième ordre est l'équation suivante.

f(x) = b_0 + b_1 x + b_2 x^2 + b_3 x^3 + b_4 x^4 + b_5 x^5

Tout d'abord, les résultats sont présentés ci-dessous. La ligne verte correspond à la fonction sin et le point bleu correspond aux données utilisées pour l'apprentissage. De plus, le rouge est le résultat de l'ajustement simple d'un polypole du cinquième ordre par la méthode des moindres carrés.

Le programme est le suivant. L'ajustement par la méthode des moindres carrés est effectué à l'aide de scipy.


import numpy as np
from scipy.odr import odrpack as odr
from scipy.odr import models
import matplotlib.pyplot as plt
import theano.tensor as tt
import scipy.stats as stat
import pymc3 as pm

#Fonction Sin
np.random.seed(5)
x = np.random.uniform(0, 10, size=20)
y = np.random.normal(np.sin(x), 0.2) #distribution normale, (moyenne,Distribué)

#Montage par scipy
polyFunc = models.polynomial(5) #5e
data = odr.Data(x, y)
myodr = odr.ODR(data, polyFunc, maxit=200)

myodr.set_job(fit_type=2) #Méthode du carré minimum

fit = myodr.run()

coeff = fit.beta[::-1]
err = fit.sd_beta[::-1]
print(coeff) #[-2.48756876e-05 -1.07573181e-02  2.09111045e-01 -1.21915241e+00, 2.11555200e+00 -1.08779827e-01] #x50

Régression polygonale par Bayes

Dans ce chapitre, sur la base des données générées dans le chapitre précédent, nous trouverons un polypole du cinquième ordre dans le cadre bayésien. Le modèle est comme le programme ci-dessous

Moyenne générée à partir d'un polypole avec des paramètres générés à partir d'une distribution normale
Distribuez les données générées par la distribution Koshi
Degré de liberté pour les données générées à partir de fonctions exponentielles (paramètres qui contrôlent la façon dont les données sont centrées dans la distribution normale)

La distribution t de l'étudiant

{\rm St}(x|\mu,\lambda,\nu)={\Gamma(\nu/2+1/2)\over\Gamma(\nu/2)}\left({\lambda\over\pi\nu}\right)^{1/2}\left[1+{\lambda(x-\mu)^2\over\nu}\right]^{-\nu/2-1/2}

Et. Ensuite, les paramètres sont appris en utilisant les méthodes de Monte Carlo en chaîne de Markov (MCMC).

with pm.Model() as model:
  b0 = pm.Normal('b0', mu=0, sd=10)
  b1 = pm.Normal('b1', mu=0, sd=5)
  b2 = pm.Normal('b2', mu=0, sd=5)
  b3 = pm.Normal('b3', mu=0, sd=5)
  b4 = pm.Normal('b4', mu=0, sd=5)
  b5 = pm.Normal('b5', mu=0, sd=5)

  #Distribution de Cauchy x ∈[0, ∞) https://docs.pymc.io/api/distributions/continuous.html#pymc3.distributions.continuous.HalfCauchy
  sd = pm.HalfCauchy('sd', 5)   
  mu = b0 + b1 * x + b2 * x**2 + b3 * x**3 + b4 * x**4 + b5 * x**5
  nu = pm.Exponential('nu', 1/30)

  y_pred = pm.StudentT('y_pred', mu=mu, sd=sd, nu=nu, observed=y)
  step = pm.Metropolis()
  trace = pm.sample(100000,step=step)

Le polypole utilisant les paramètres estimés est indiqué en rouge ci-dessous. La ligne verte est la fonction sin, le point bleu est les données utilisées pour l'entraînement et le bleu est le résultat de la méthode du carré minimum. De plus, la moyenne des variances obtenues à la suite du MCMC était de 0,23, et la variance du modèle généré pouvait être estimée. Cependant, on peut voir que l'erreur est grande à partir de la fonction sin réelle, et même dans la partie où il y a peu de données, comme autour de x = 10, la variance devient un polynôme généré à partir d'une constante 0,23. C'est une mauvaise conception du modèle.

Développement dans le processus gaussien

Dans le chapitre précédent, la fonction f (x) était représentée par un modèle polymorphe et ses paramètres ont été déduits. Cependant, les données ne sont pas bien représentées et il y a de fortes chances que le modèle soit erroné. Par conséquent, nous utilisons la méthode pour déduire directement la fonction. Autrement dit, dans le chapitre précédent, p(θ|x)L'estimation des paramètres a été effectuée, mais dans ce chapitre, p(f|x)Déduire. Processus gaussien comme méthode pour déduire cette fonction(GP)Utiliser.

La définition officielle de GP est la suivante. Tous les points dans l'espace continu sont associés à des variables qui suivent une distribution normale, et GP est la distribution simultanée de leur nombre infiniment grand de variables stochastiques.

Mais je ne vais pas avoir une discussion particulièrement difficile ici. Il y a trois choses importantes:

La distribution a priori est une distribution normale multivariée. f (x) ~ MvNormal (u | K)
La matrice de covariance (K) est le noyau.
Paramètres de la matrice de covariance
La distribution postérieure peut être calculée analytiquement.

Imaginez le noyau de 2 utilisé dans une machine à vecteurs de support ou autre. Il existe différents noyaux et n'importe quel noyau peut être utilisé comme matrice de codistribution. Cet article utilise le noyau gaussien. Jetons un coup d'œil à la distribution antérieure de GP obtenue à partir de 1 et 2. Le programme suivant trace une distribution normale multivariée (moyenne = 0, matrice de covariance = noyau gaussien). En regardant la figure, vous pouvez voir qu'elle représente diverses fonctions. Il existe d'innombrables expressions comme celle-ci.

Notez que si vous utilisez un noyau différent, les expressions telles que les fonctions irrégulières changeront. Cependant, cette fois, lorsqu'on considère la régression de la fonction sin, on peut voir que l'expression montrée dans la figure ci-dessous est appropriée.

def squareDistance(x, y):
  return np.array([[(x[i]-y[j])**2 for i in range(len(x))] for j in range(len(y))])

np.random.seed(1)
tp = np.linspace(0, 10, 100)

cov = np.exp(-squareDistance(tp, tp))
plt.plot(tp, stats.multivariate_normal.rvs(cov=cov, size=6).T

La distribution normale multivariée illustrée dans la figure ci-dessus est exprimée de diverses manières en modifiant la matrice de covariance. Afin d'estimer la matrice de covariance dans le cadre bayésien, il est nécessaire de définir les paramètres à estimer comme indiqué en 3. Donc, la matrice de covariance

K_{i,j} = \eta \exp(-\rho D)  \ \  (i \neq j) \\
K_{i,j} = \eta + \sigma \ \ (i = j)

ça ira.

Enfin, il existe quatre calculs de distribution a posteriori analytiques. Si la probabilité est une distribution normale, la distribution postérieure GP peut résoudre le problème de manière analytique, de sorte que la distribution normale multivariée, qui est la distribution postérieure GP, peut être calculée avec une formule simple. La distribution postérieure permet de prédire la fonction de régression.

Une explication détaillée de la formule peut être trouvée dans [Gaussian Process and Machine Learning] Gaussian Process Return.

Les résultats de la prédiction et une partie du programme sont présentés ci-dessous. Le résultat prédit à partir de la distribution postérieure estimée de GP est la ligne rouge, la ligne verte est la fonction sin, le point bleu est les données utilisées pour l'entraînement et la ligne bleue est le résultat de la méthode du carré minimum. En raison de la prédiction par GP, la ligne rouge est plus proche de la ligne verte que de la ligne bleue, on peut donc voir qu'un modèle de régression avec de bonnes performances peut être estimé directement sans ajuster les paramètres. En outre, en regardant comment les lignes rouges sont dispersées, on peut voir que les lignes sont dispersées, c'est-à-dire que la dispersion est grande, là où il n'y a pas de données telles que 10 environ.

#Distribution antérieure GP: f(x) ~ GP(u=(u1, u2, ...), k(x, x')), u1, u2 .. =Supposé 0
#Responsabilité: p(y|x, f(x)) ~ N(f, σ^2I)
#Distribution postérieure GP p(f(x)|x,y) ~GP(U,∑)

def squareDistance(x, y):
  return np.array([[(x[i]-y[j])**2 for i in range(len(x))] for j in range(len(y))])

points = np.linspace(xmin, xmax, xsize) #Point de référence

with pm.Model() as model:
  #Distribution antérieure
  mu = np.zeros(xsize)
  eta = pm.HalfCauchy('eta', 5)
  rho = pm.HalfCauchy('rho', 5)
  sigma = pm.HalfCauchy('sigma', 5)

  D = squareDistance(x, x)
  K = tt.fill_diagonal(eta * pm.math.exp(-rho * D), eta + sigma) #i=eta en position j+Stocker sigma

  gpPri = pm.MvNormal('obs', mu, cov=K, observed=y)

  #Calcul de la distribution postérieure(La distribution postérieure de la distribution gaussienne peut être calculée par la formule)
  K_ss = eta * pm.math.exp(-rho * squareDistance(points, points))
  K_s = eta * pm.math.exp(-rho * squareDistance(x, points))

  MuPost = pm.Deterministic('muPost', pm.math.dot(pm.math.dot(K_s, tt.nlinalg.matrix_inverse(K)), y)) #K_s.T inv(K) y
  SigmaPost = pm.Deterministic('sigmaPost', K_ss - pm.math.dot(pm.math.dot(K_s, tt.nlinalg.matrix_inverse(K)), K_s.T))

  step = pm.Metropolis()
  trace = pm.sample(10000,step=step)

Résumé

Dans cet article, nous avons présenté les statistiques bayésiennes car l'expression de l'incertitude est importante. Cependant, lors de la création d'un modèle de régression de la fonction sin, il était difficile de concevoir un modèle qui exprime bien la fonction sin, nous avons donc introduit le processus gaussien. J'éviterai les explications détaillées et expliquerai le concept. Pour plus de détails, veuillez consulter les autres articles et livres de l'article.

[PYTHON] Retour en utilisant le processus gaussien

Aperçu

Régression polymérique de la fonction Sin

Régression polygonale par Bayes

Développement dans le processus gaussien

Résumé