Processus gaussien avec pymc3

J'essaierai d'estimer par Gaussian Process en utilisant les données sur les ventes de shampooing trouvées sur ce site.

https://machinelearningmastery.com/time-series-datasets-for-machine-learning/

Le processus gaussien est un modèle de régression non paramétrique typique. Il présente l'avantage d'être solide par rapport aux données non linéaires et de pouvoir confirmer l'incertitude de l'estimation.

import os
import pandas as pd
shampoo_df = pd.read_csv('../data/shampoo.csv')
shampoo_df
shampoo_df.plot.line(x='Month', y='Sales',c="k");

Dans ce qui suit, la vraisemblance marginale est calculée afin d'obtenir une vraisemblance appropriée et une paire de distribution antérieure.

with pm.Model() as shampoo_model:
    ρ=pm.HalfCauchy('ρ', 1) #Plus l'échelle de longueur est grande, plus les mouvements semblent proches.<= Auto-corrélé
    η = pm.HalfCauchy('η', 1) #Plus l'échelle de longueur est grande, plus les mouvements semblent proches.<= Auto-corrélé
    
    M = pm.gp.mean.Linear(coeffs=(shampoo_df.index/shampoo_df.Sales).mean()) #valeur initiale
    K= (η**2) * pm.gp.cov.ExpQuad(1, ρ)

    σ = pm.HalfCauchy('σ', 1)
    gp = pm.gp.Marginal(mean_func=M, cov_func=K) 
    gp.marginal_likelihood("recruits", X=shampoo_df.index.values.reshape(-1,1), 
                           y=shampoo_df.Sales.values, noise=σ)
    trace = pm.sample(1000)

pm.traceplot(trace);

L'estimation la plus probable crée un modèle bien ajusté en maximisant la fonction de vraisemblance p (X | θ) par rapport au paramètre θ, mais n'utilise pas la distribution antérieure p (θ). D'un autre côté, il semble que le cas où il y a une distribution a priori s'appelle l'estimation de probabilité postérieure maximale.

Veuillez vous référer aux livres suivants pour plus de détails. Atsushi Suyama. Série de démarrage d'apprentissage automatique Introduction à l'apprentissage automatique par l'inférence bayésienne (édition japonaise) (emplacements Kindle 2154-2156). Édition Kindle.

Obtenez les paramètres optimaux en estimant la probabilité postérieure maximale.

with shampoo_model:
    marginal_post = pm.find_MAP()

{'ρ_log__': array(2.88298779), 'η_log__': array(5.63830367), 'σ_log__': array(4.13876033), 'ρ': array(17.86757799), 'η': array(280.98567051), 'σ': array(62.72501496)}

X_pred=np.array([r for r in range(X_pred.max()+10)]).reshape(-1,1)

with shampoo_model:
    shampoo_pred = gp.conditional("shampoo_", X_pred)
    samples=pm.sample_ppc([marginal_post], vars=[shampoo_pred] ,samples=100)



fig = plt.figure(figsize=(10,5)); 
ax = fig.gca()

from pymc3.gp.util import plot_gp_dist
plot_gp_dist(ax, samples["shampoo_"], X_pred);

ax.plot(shampoo_df.index, shampoo_df.Sales, 'dodgerblue', ms=5, alpha=0.5);
ax.plot(shampoo_df.index, shampoo_df.Sales, 'ok', ms=5, alpha=0.5);
plt.show()

Dans pm.sample_ppc, l'échantillonnage est effectué à partir de la distribution de prédiction postérieure basée sur les paramètres de marginal_post.