introduction

J'avais peur qu'il n'y ait pas beaucoup d'implémentations correctes de la régression linéaire bayésienne dans le monde, et il y a peu d'implémentations qui prennent en charge l'entrée multidimensionnelle, donc je l'ai implémentée comme une classe facile à utiliser. La description suit essentiellement PRML.

Classe implémentée

Cela dit, ce n'est pas aussi gros qu'une implémentation, c'est environ 50 lignes ... J'écrirai toute la classe ici. Ses fonctions sont les suivantes. C'est très simple.

-Donner une combinaison de $ \ phi $ et $ t $ pour calculer la distribution postérieure

Échantillonnage aléatoire des paramètres de moyenne et de covariance pour la régression linéaire bayésienne -Donner $ \ phi $ pour sortir la moyenne et la variance de la distribution prédite

`Beyes_LR.py`


import numpy as np
from scipy.stats import multivariate_normal

class BeyesLinearRegression:
    def __init__(self, mu, S, beta):
        self.mu = mu
        self.S = S
        self.beta = beta

    def calc_posterior(self, phi, t):
        S_inv = np.linalg.inv(self.S)

        if len(phi.shape) == 1:
            phi = phi.reshape(1, -1)
            t = t.reshape(1, 1)
        self.S = np.linalg.inv(S_inv + self.beta * phi.T @ phi)
        self.mu = self.S @ (S_inv @ self.mu + np.squeeze(self.beta * phi.T @ t))

    def sampling_params(self, n=1, random_state=0):
        np.random.seed(random_state)
        return np.random.multivariate_normal(self.mu, self.S, n)

    def probability(self, x):
        dist = multivariate_normal(mean=self.mu, cov=self.S)
        return dist.logpdf(x)

    def predict(self, phi):
        if len(phi.shape) == 1:
            phi = phi.reshape(1, -1)
        pred = np.array([self.mu.T @ _phi for _phi in phi])
        S_pred = np.array([(1 / self.beta) + _phi.T @ self.S @ _phi for _phi in phi])

        # Above is a simple implementation.
        # This may be better if you want speed.
        # pred = self.mu @ phi.T
        # S_pred = (1 / self.beta) + np.diag(phi @ self.S @ phi.T)
        return pred, S_pred

Tout le code est également dans git. (Bien que ce soit environ 50 lignes)

GitHub

À propos de la régression linéaire bayésienne

Puisque ce qui suit est détaillé sur la dérivation de la formule, je n'écrirai pas les détails.

Calcul de la distribution postérieure

L'important est de mettre à jour la distribution comme suit. Pour le dire clairement, $ \ phi $ est la variable explicative et $ t $ est la réponse. Les matrices de moyenne et de covariance sont mises à jour en conséquence.

M_N=S_N(S_0^{-1}m_0+\beta\Phi^Tt)

S_N^{-1}=S_0^{-1}+\beta\Phi^T\Phi

Distribution prévue

La distribution prévue est indiquée ci-dessous. Je vais également omettre les détails ici, mais le fait est que la distribution est exprimée par la moyenne et la variance comme indiqué ci-dessous pour le nouveau point $ x $.

N(m_N^T\phi(x), 1/\beta+\phi(x)^TS_N\phi(x))

Utilisez la classe implémentée

À partir de là, essayons la régression linéaire bayésienne en utilisant la classe que nous avons réellement implémentée. La classe que j'ai créée doit recevoir $ \ phi $ en tant que fonctionnalité de $ x $. Dans une implémentation courante, la partie génération de $ \ phi $ (par exemple, un polynôme) est également incluse dans la classe, et il est difficile de dire s'il exécute une régression linéaire bayésienne ou s'il conçoit des entités à l'aide de polynômes. Mais ici, il est séparé.

Donc, si vous entrez les données originales $ x $, $ y $ puis implémentez la fonction pour créer $ \ phi $, c'est fondamentalement OK.

Régression avec données d'onde sinusoïdale

Essayons d'abord avec les données Toy.

Génération de données et conception de fonctions

Les données d'entrée sont les données d'observation en ajoutant du bruit à la distribution réelle de l'onde sinusoïdale. De plus, la fonction caractéristique est conçue comme une onde composite de fonctions triangulaires de plusieurs fréquences. La méthode x_to_phi vectorise 10 ondes et _phi représente une onde composite. L'amplitude est le paramètre obtenu par régression linéaire bayésienne. L'image ci-dessous est mathématique. (Si vous y réfléchissez, le premier élément est zéro ... je n'en ai pas besoin ...)

y=w_1sin(0)+w_2sin(2\pi x)+w_3sin(2\times2\pi x)+\cdots+w_9sin(9\times2\pi x)

def x_to_phi(x):
    if len(x.shape) == 1:
        x = x.reshape(-1, 1)
    return np.concatenate([np.sin(2 * np.pi * x * m) for m in range(0, 10)], axis=1)


def _phi(x, params):
    return np.array([p * np.sin(2 * np.pi * x * i) for i, p in enumerate(params)]).sum(axis=0)

Cliquez ici pour la partie réelle de la génération de données.

x = np.arange(0, 1, 0.01)
phi = x_to_phi(x)

e = np.random.randn(len(x)) / 10
y_true = np.sin(2 * np.pi * x)
y = y_true + e

Lorsqu'un seul point est observé

Tout d'abord, considérons le cas où un seul point des 50e données est observé.

train_idx = 50
x_train = x[train_idx]
phi_train = phi[train_idx]
y_train = y[train_idx]
plt.scatter(x_train, y_train, c='crimson', marker='o', label='observation')
plt.plot(x, y_true, label='true')

Calculons immédiatement la distribution a posteriori de ces données. Si vous voulez juste apprendre, c'est une ligne comme celle-ci:

#Valeur initiale de la régression linéaire bayésienne
input_dim = phi.shape[1]
mu = np.zeros(input_dim)
S = np.identity(input_dim)
sigma = 0.1
beta = 1.0 / (sigma ** 2)

#Définition du modèle et formation
beyes_linear_model = BeyesLR.BeyesLinearRegression(mu, S, beta)
beyes_linear_model.calc_posterior(phi_train, y_train)

Certaines formes d'ondes échantillonnées au hasard de la distribution postérieure post-entraînement sont affichées avec des lignes pointillées vertes et la distribution prédite est affichée en bleu clair. La plupart d'entre eux sont bleus car je n'ai appris qu'un seul point.

Lors de l'observation de 50 points

Faisons exactement la même chose à partir de 50 données d'observation. La seule différence dans le code est de sélectionner au hasard 50 train_idx. La distribution prévue des résultats est indiquée dans la figure ci-dessous.

Retour dans l'ensemble de données publicitaire

Vient ensuite un vrai problème, qui résout également une régression linéaire multidimensionnelle. Si vous extrayez des entités dans plusieurs dimensions, les dimensions augmenteront trop et cela sera difficile à comprendre, donc $ \ phi $ est une fonction de Linear.

Des données d'entrée

Il s'agit d'un ensemble de données publicitaires familier à ISLR. Cette fois, nous utiliserons les dépenses publicitaires TV et Radio comme intrants et les ventes comme réponse.

Cette fois, $ \ phi $ est linéaire, alors ajoutez simplement le terme de section. La formule de la régression est la suivante. $Sales = w_0+w_1TV+w_2Radio$

def x_to_phi(x, typ='linear', degree=3):
    if len(x.shape) == 1:
        x = x.reshape(-1, 1)
    return np.concatenate([np.ones(x.shape[0]).reshape(-1, 1), x], axis=1)


df = pd.read_csv(ADVERTISING_DATASET)
x = df[['TV', 'Radio']].values
y = df['Sales'].values

phi = x_to_phi(x)
x_train, x_test, phi_train, phi_test, y_train, y_test = \
    train_test_split(x, phi, y, train_size=0.05, random_state=0)

Tout ce que vous avez à faire est d'apprendre comme dans l'exemple précédent.

input_dim = phi.shape[1]
mu = np.zeros(input_dim)
S = np.identity(input_dim)
sigma = 10
beta = 1.0 / (sigma ** 2)

beyes_linear_model = BeyesLR.BeyesLinearRegression(mu, S, beta)
beyes_linear_model.calc_posterior(phi_train, y_train)

Dans le code, train_size est défini sur 0,05, mais le plan de régression dessiné lors de sa modification est le suivant. Il est appris par régression linéaire bayésienne, et 5 plans sont extraits et dessinés par échantillonnage aléatoire. Aléatoire lorsque le nombre d'échantillons d'apprentissage est petit Un plan est dessiné, mais il converge à mesure que le nombre de points de données augmente.

en conclusion

Enfin, une petite promotion de la régression linéaire bayésienne. Bien que certaines parties ne soient toujours pas entièrement comprises, la conception de la fonction d'extraction de caractéristiques $ \ phi $ de la régression linéaire bayésienne devient importante. Je reconnais que la régression de processus gaussienne traite la matrice distribuée co-distribuée comme une matrice de planification utilisant des fonctions du noyau sans l'écrire explicitement. Cependant, l'expérience montre que la régression linéaire bayésienne est suffisante du point de vue de la descriptivité s'il existe des connaissances préalables capables d'extraire de bonnes caractéristiques. De plus, grâce à la méthode de calcul de la régression linéaire bayésienne, l'apprentissage séquentiel peut être effectué tel quel. Tout ce que vous avez à faire est d'apprendre la distribution postérieure calculée en tant que distribution antérieure. Il n'est pas nécessaire de recalculer la matrice de gramme comme la régression de processus gaussien. Il existe peut-être une version en ligne, mais ...

Au fait, ayez une vie de retour linéaire bayésienne confortable!

Implémentation python de la classe de régression linéaire bayésienne