Technologie

Régression bayésienne avec implémentation dans R

Régression bayésienne avec implémentation dans R


Dérivations théoriques à partir de zéro, implémentation R et discussion de la vision du monde bayésienne

2*SOJ0 TviGc gPnHrbeOY3w - Régression bayésienne avec implémentation dans R
Un modèle graphique probabiliste montrant les dépendances entre les variables en régression (Bishop 2006)

La régression linéaire peut être établie et interprétée dans une perspective bayésienne. En écrivant sur la régression linéaire bayésienne, qui est elle-même intéressante à penser, je peux également discuter de la vision générale du monde bayésien. Les premières parties discutent de la théorie et des hypothèses à partir de zéro, et les parties suivantes incluent une implémentation R et des remarques. Les lecteurs peuvent se sentir libres de copier les deux blocs de code dans un carnet R et de jouer avec. Les commentaires sur tout ce qui est discuté ici, en particulier la philosophie bayésienne, sont plus que bienvenus.

Rappelons qu'en régression linéaire, on nous donne des valeurs cibles y, Les données X, et nous utilisons le modèle

où y est le vecteur N * 1, X est la matrice N * D, w est le vecteur D * 1 et l'erreur est le vecteur N * 1. Nous avons N points de données. La dimension D est comprise en termes de caractéristiques, donc si nous utilisons une liste de x, une liste de x² (et une liste de 1 correspondant à w_0), nous disons D = 3. Si vous n'aimez pas la forme matricielle, pensez-y simplement comme une forme condensée des éléments suivants, où tout est un nombre au lieu d'un vecteur ou d'une matrice:

Dans la régression linéaire classique, le terme d'erreur est supposé avoir une distribution normale, et il s'ensuit immédiatement que y est normalement distribué avec une moyenne Xw, et la variance de la variance du terme d'erreur (dénotée par σ² ou matrice diagonale avec les entrées σ²). L'hypothèse normale s'avère bien dans la plupart des cas, et ce modèle normal est également ce que nous utilisons dans la régression bayésienne.

Nous sommes maintenant confrontés à deux problèmes: l'inférence de w et la prédiction de y pour tout nouveau X. En utilisant la règle bien connue de Bayes et les hypothèses ci-dessus, nous ne sommes qu'à quelques pas non seulement de résoudre ces deux problèmes, mais distribution de probabilité complète de y pour tout nouveau X. Voici la règle de Bayes en utilisant nos notations, qui exprime la distribution postérieure du paramètre w données données:

π et f sont des fonctions de densité de probabilité. Puisque le résultat est une fonction de w, nous pouvons ignorer le dénominateur, sachant que le numérateur est proportionnel au côté gauche par une constante. Nous savons par hypothèses que la fonction de vraisemblance f (y | w, x) suit la distribution normale. L'autre terme est la distribution préalable de w, et cela reflète, comme son nom l'indique, une connaissance préalable des paramètres.

Distribution préalable Définir le prieur est une partie intéressante du flux de travail bayésien. Pour plus de commodité, nous laissons w ~ N (m_0, S_0), et les hyperparamètres m et S reflètent maintenant la connaissance antérieure de w. Si vous avez peu de connaissances sur w, ou trouvez une affectation de m et S trop subjective, les prieurs «non informatifs» sont un amendement. Dans ce cas, nous mettons m à 0 et plus important encore, S comme matrice diagonale avec de très grandes valeurs. Nous disons que w a une variance très élevée, et nous avons donc peu de connaissances sur ce que sera w.

Avec toutes ces fonctions de probabilité définies, quelques lignes de manipulations simplement algébriques (pas mal de lignes en fait) donneront le postérieur après observation de N points de données:

Cela ressemble à un tas de symboles, mais ils sont déjà tous définis, et vous pouvez calculer cette distribution une fois que ce résultat théorique est implémenté dans le code. (N (m, S) signifie une distribution normale avec une moyenne m et une matrice de covariance S.)

Une approche bayésienne complète signifie non seulement obtenir une seule prédiction (dénoter une nouvelle paire de données par y_0, x_0), mais également acquérir la distribution de ce nouveau point.

Ce que nous avons fait est l'inverse de la marginalisation du joint pour obtenir la distribution marginale sur la première ligne, et en utilisant la règle de Bayes à l'intérieur de l'intégrale sur la deuxième ligne, où nous avons également supprimé les dépendances inutiles. Notez que nous savons quelles sont les deux dernières fonctions de probabilité. Le résultat d'une distribution prédictive complète est:

L'implémentation en R est assez pratique. Soutenu par les résultats théoriques ci-dessus, nous entrons simplement des multiplications matricielles dans notre code et obtenons les résultats des prédictions et des distributions prédictives. Pour illustrer avec un exemple, nous utilisons un problème de jouet: X est de -1 à 1, régulièrement espacé, et y est construit comme les ajouts suivants de courbes sinusoïdales avec un bruit normal (voir les images ci-dessous pour une illustration de y).

Le code suivant obtient ces données.

bibliothèque (ggplot2)

# - - - - - Obtenir des données - - - - - - - - - - - - - - - - - - - - - +

X <- (-30: 30) / 30
N <- longueur (X)
D <- 10
var <- 0,15 * 0,15
e <- rnorm (N, 0, var ^ 0,5)
EY <- sin (2 * pi * X) * (X0)
Y <- sin (2 * pi * X) * (X0) + e
data <- data.frame (X, Y)
g1 <- ggplot (data = data) + geom_point (mapping = aes (x = X, y = Y))

Le code suivant (sous la section Maximum a Posteriori) implémente tous les résultats théoriques ci-dessus par multiplications matricielles. Nous développons également les fonctionnalités de x (désignées dans le code comme phi_X, dans la section Construire les fonctions de base). Tout comme nous étendrions x à x², etc., nous l'étendons maintenant en 9 fonctions de base radiales, chacune ressemblant à ce qui suit. Notez que bien qu'elles ressemblent à une densité normale, elles ne sont pas interprétées comme des probabilités.

Un avantage des fonctions de base radiales est que les fonctions de base radiales peuvent s'adapter à une variété de courbes, y compris polynomiales et sinusoïdales.

# - - - - - Construire des fonctions de base - - - - - - - - - - - - +phi_X <- matrice (0, nrow = N, ncol = D)
phi_X[,1] <- X
mu <- seq (min (X), max (X), longueur.out = D + 1)
mu <- mu[c(-1,-length(mu))]
pour (i dans 2: D) {
phi_X[,i] <- exp (- (X-mu[i-1]) ^ 2 / (2 * var))
}
# - - - - - Maximum a Posteriori - - - - - - - - - - - - - - - +# Les commentaires sont généraux avant
# m0 <- matrice (0, D, 1)
# S0 <- diag (x = 1000, D, D)
# SN <- inv (inv (S0) + t (phi_X)% *% phi_X / var)
# mN <- SN% *% (inv (S0)% *% m0 + t (phi_X)% *% Y / var)
# Y_hat <- t (mN)% *% t (phi_X)
# Nous utilisons avant non informatif pour l'instant
m0 <- matrice (0, D, 1)
SN <- résoudre (t (phi_X)% *% phi_X / var)
mN <- SN% *% t (phi_X)% *% Y / var
Y_hat <- t (mN)% *% t (phi_X)
var_hat <- tableau (0, N)
pour (i en 1: N) {
var_hat[i] <- var + phi_X[i,]% *% SN% *% phi_X[i,]
}
g_bayes <- g1 +
geom_line (mapping = aes (x = X, y = Y_hat[1,]), color = ’# 0000FF’)
g_bayes_full <- g_bayes + geom_ribbon (mapping = aes (x = X, y = Y_hat[1,],
ymin = Y_hat[1,]-1,96 * var_hat ^ 0,5,
ymax = Y_hat[1,]+ 1,96 * var_hat ^ 0,5, alpha = 0,1),
fill = ’# 9999FF’)

Un détail à noter dans ces calculs, c'est que nous utilisons au préalable des informations non informatives. La section commentée est exactement les résultats théoriques ci-dessus, tandis que pour les informations non informatives antérieures, nous utilisons une matrice de covariance avec des entrées diagonales approchant l'infini, donc l'inverse de celui-ci est directement considéré comme 0 dans ce code. Si vous souhaitez utiliser ce code, assurez-vous d'installer le package ggplot2 pour le traçage.

L'illustration suivante vise à représenter une distribution prédictive complète et à donner une idée de l'adéquation des données.

La ligne bleue est la valeur attendue de la distribution prédictive à chaque point x, et la région bleu clair se réfère aux régions à l'intérieur de deux écarts-types. La ligne rouge est la véritable fonction de y. Les points sont des données générées aléatoirement à partir de la fonction donnée avec un bruit normal.

Le résultat de la régression linéaire multiple est identique à celui de la régression bayésienne utilisant un a priori impropre avec une matrice de covariance infinie. Cependant, la distribution prédictive de la régression bayésienne présente généralement une variance plus étroite. En règle générale, il est recommandé d'obtenir une certaine connaissance du domaine concernant les paramètres et d'utiliser un préalable informatif. La régression bayésienne peut ensuite rapidement quantifier et montrer comment les différentes connaissances antérieures ont un impact sur les prédictions.

La régression bayésienne est assez flexible car elle quantifie toutes les incertitudes - les prévisions et tous les paramètres. Cette flexibilité offre plusieurs commodités. Par exemple, vous pouvez marginaliser toutes les variables des distributions conjointes et étudier la distribution de toutes les combinaisons de variables. De plus, l'ajustement des données dans cette perspective vous permet «d'apprendre au fur et à mesure». Disons que j'ai d'abord observé 10000 points de données et calculé une valeur postérieure du paramètre w. Après cela, j'ai en quelque sorte réussi à acquérir 1000 points de données supplémentaires, et au lieu d'exécuter à nouveau toute la régression, je peux utiliser le postérieur précédemment calculé comme mon précédent pour ces 1000 points. Ce processus séquentiel donne le même résultat que d'utiliser à nouveau toutes les données. J'aime cette idée dans la mesure où elle est très intuitive, dans la mesure où une opinion apprise est proportionnelle à des opinions précédemment apprises et à de nouvelles observations, et l'apprentissage se poursuit. Une plaisanterie dit qu'un Bayésien qui rêve d'un cheval et observe un âne, l'appellera un mulet. Mais s'il en prend plus, il finira par dire qu'il s'agit bien d'un âne.

Afficher plus

SupportIvy

SupportIvy.com : Un lieu pour partager le savoir et mieux comprendre le monde. Meilleure plate-forme de support gratuit pour vous, Documentation &Tutoriels par les experts.

Articles similaires

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *

Bouton retour en haut de la page
Fermer