[PYTHON] Essayez le modèle d'espace d'état (Jupyter Notebook + noyau IR)

Aperçu du modèle d'espace d'états

Le livre est Iwanami Data Science Vol.1. ([Dossier spécial] Raisonnement bayésien et logiciel gratuit MCMC) Le diagramme conceptuel du modèle d'espace d'états est cité ici.

Le modèle State-Space a une série chronologique ** «état» **, et sur cette base, une ** «valeur observée» ** contenant une certaine variation est générée. C'est. Iwanami Data Science Vol.1 est une fonction spéciale de "Bayes Inference and MCMC Free Software", mais le chapitre 3 est un article sur les modèles d'espace d'état (écrit par M. Ito), et d'autres chapitres expliquent également les modèles spatiaux. viens. Le chapitre 3 a expliqué comment résoudre ce problème avec R et JAGS (Just Another Gibbs Sampler).

Quand j'ai essayé cela, je l'ai essayé avec Python + PyMC3 au début, mais cela n'a pas fonctionné.

import pymc3
(Omis)
model1 = pm.Model()

with model1:  # model definition
    alpha0 = pm.Uniform('alpha0', lower=0. , upper=10.)
    sigma_y = pm.Uniform('sigma_y', lower=0., upper=100.)
    sigma_alp = pm.Uniform('sigma_alp', lower=0., upper=100.)
    
    for i in range(1,N):   #Les deux lignes suivantes sont ERROR.Je voulais définir la quantité d'état α...
        alpha[i] = pm.Normal(1, mu=alpha[i-1], sd=sigma[0])

    y_obs = pm.Normal('y_obs', mu=alpha, sd=sigma[0], observed=ydat)
(Omis)

J'ai essayé ce genre de code, mais dans la boucle for, j'obtiens une erreur au point où le modèle de quantité d'état $ \ alpha_t = \ alpha_ {t-1} + \ eta $ est défini. Si vous y réfléchissez bien, ** pymc3.Normal () ** (** pm.Normal () ** dans la liste) est un constructeur de classe pour générer un modèle de distribution normale, donc dans son argument Il est considéré difficile d'inclure la classe (alpha [i-1]) à générer. (J'ai cherché sur Internet de diverses manières, mais je n'ai pas trouvé de solution.)

J'ai envisagé d'utiliser ** PyStan ** et l'implémentation du système de filtrage de Kalman ** PyKalman ** comme d'autres bibliothèques capables de gérer de tels modèles d'espace d'état (ou modèles linéaires dynamiques), mais cette fois j'ai suivi le livre J'ai décidé d'utiliser "R". Le programme utilisé est le suivant.

• R ver.3.2.2 --Jupyter Notebook (notebook est ver.4.0.6) + IRkernel ... Cette fois, j'ai essayé IRkernel pour la première fois. --JAGS ver.4.0.0 + {rjags} package
• Paquet {dlm}

À propos de l'environnement de Jupyter Notebook + IRkernel

Jupyter Notebook est une version mise à jour du notebook IPython que vous connaissez. En commutant le système de traitement du langage et le noyau, il est possible de manipuler autre chose que Python. (Il a été nommé d'après Jupyter = Julia + Python + R.) Cette fois, j'ai installé et utilisé IRkernel pour le langage R.

La procédure d'installation est la suivante.

-Installer le langage R (si nécessaire). --Installez IRkernel (https://github.com/IRkernel/IRkernel) sous l'environnement R.

install.packages(c('rzmq','repr','IRkernel','IRdisplay'),
                 repos = c('http://irkernel.github.io/', getOption('repos')))
IRkernel::installspec()

• Commencez par "> Jupyte Notebook" sur la console. (Au début, une erreur s'est produite en raison de l'absence de certains packages dépendants de R, mais elle peut être démarrée en installant le package R comme demandé par le message d'erreur.)

** Fig. Au démarrage de Jupyter Notebook **

Comme indiqué ci-dessus, vous pouvez sélectionner "R" dans le nouveau menu de création de note.

Implémentation du modèle d'espace d'états par JAGS

Suivez le chapitre «Modélisation de séries temporelles et de données spatiales» dans le livre «Iwanami Data Science Vol.1» pour créer les notes suivantes.

Comme avec les blocs-notes IPython précédents, vous pouvez mélanger du texte Markdown avec du code exécutable. De plus, à partir de cette mise à jour de Notebook, je suis heureux de pouvoir rendre les formules ** MathJax ** hors ligne.

Une astuce pour tracer une figure avec R. (Notebook In [1]: il s'agit de la commande décrite.)

library(repr)

# Change plot size to 4 x 3
options(repr.plot.width=4, repr.plot.height=3)

Dans la seconde moitié de cette note, nous traçons en utilisant {ggplot2}, mais si rien n'est spécifié, la figure sera générée pour s'adapter à toute la largeur du navigateur Web. Le chiffre est assez grand (bien que la plupart des gens se trouvent dans une situation similaire), et l'équilibre n'est pas très bon lorsqu'il est affiché avec des phrases et des codes. {repr} est un sous-paquet qui est installé avec IRkernerl, mais vous pouvez utiliser cette fonction pour spécifier des options comme ci-dessus pour créer la figure de n'importe quelle taille.

De plus, bien que cette note n'exécute pas directement le code, le code BUGS est également décrit. (JAGS + {rjags} est une spécification qui prépare le code BUGS dans un fichier séparé et exécute le calcul.)

J'ai senti que c'était un grand avantage d'utiliser Jupyter Notebook pour pouvoir combiner le code R qui résume le tout, le code BUGS qui définit le modèle de simulation, le résultat du calcul et les informations qui ont tendance à être dispersées séparément.

La transition de la largeur annuelle de l'anneau de Sugi à Kyoto calculée par le modèle de niveau local de cette manière est illustrée comme suit.

Fig. 'NenRin' growth width (by JAGS)

La ligne pointillée avec le marqueur est la valeur observée, la ligne pleine épaisse est la valeur estimée par le calcul de MCMC, et la zone grise est la section de crédit de 95%. (Le code R peut être obtenu sur le site d'assistance "Iwanami Data Science".)

Confirmation par le package {dlm}

JAGS est un système de mise en œuvre de MCMC (chaîne Marcov Monte Carlo). En plus de MCMC, il existe une implémentation du système de filtre de Kalman comme méthode de calcul qui peut gérer les modèles d'espace d'états (modèles linéaires dynamiques), mais dans "R", les packages {dlm} et {KFAS} sont souvent utilisés. Cette fois, j'ai essayé le package {dlm}.

La différence entre MCMC et filtre de Kalman est expliquée correctement dans "Iwanami Data Science". (Je vais citer ...) "La différence entre le filtre de Kalman et MCMC est que la distribution postérieure périphérique peut être calculée pour l'état $ x_i $, mais elle ressemble à $ s ^ 2 $, $ \ sigma ^ 2 $. Pour les paramètres du modèle, l'estimation (estimation ponctuelle) est effectuée en ignorant l'étalement de la distribution postérieure. »L'étalement des paramètres de la relation de dispersion, pour ainsi dire, la« variation de variation »n'est pas très intéressé (personnellement). ， Ce point n'a pas d'importance au début.

Ci-dessous, le code R est cité à partir du bloc-notes créé.

## Year and Growth width data(mm)
X <- 1961:1990
Y <- c(4.71, 7.70, 7.97, 8.35, 5.70,
       7.33, 3.10, 4.98, 3.75, 3.35,
       1.84, 3.28, 2.77, 2.72, 2.54,
       3.23, 2.45, 1.90, 2.56, 2.12,
       1.78, 3.18, 2.64, 1.86, 1.69,
       0.81, 1.02, 1.40, 1.31, 1.57)

tsData <- ts(Y, frequency=1, start=c(1961, 1))   # Time Series data structure
n <- length(tsData)

library(dlm)

# Modeling function and initial parameters
mkModel1 <- function(parm){
    dlmModPoly(order = 1, dV=parm[1], dW=parm[2])
}
myparm1 <- c(0.001, 0.001)
# Parameter calculation by fitting
dlmFit1 <- dlmMLE(tsData, 
                 parm = myparm1, 
                 build = mkModel1
                 )
# Filtering
oFitted1 <- mkModel1(dlmFit1$par)
oFiltered1 <- dlmFilter(tsData, oFitted1)
print(names(oFiltered1))
print(oFiltered1$f)
# Smoothing
oSmoothed1 <- dlmSmooth(oFiltered1)
print(names(oSmoothed1))
print(oSmoothed1$s)

Avec ce qui précède, le résultat du calcul du filtre de Kalman peut être obtenu. Tracez ce résultat avec {ggplot2}.

library(ggplot2)
base.family <- ""
ribbon.alpha <- 0.3
ribbon.fill <- "black"

alpha.mean <- dropFirst(oSmoothed1$s)
var_p <- unlist(dlmSvd2var(oSmoothed1$U.S, oSmoothed1$D.S))
alpha.up <- dropFirst(oSmoothed1$s) + qnorm(0.025, sd = sqrt(var_p[-1]))
alpha.low <- dropFirst(oSmoothed1$s) + qnorm(0.975, sd = sqrt(var_p[-1]))

p3 <- ggplot(data.frame(Year = X, Width = Y,
                        Width.mean = alpha.mean,
                        Width.up = alpha.up,
                        Width.low = alpha.low)) +
  geom_ribbon(aes(x = Year, ymax = Width.up, ymin = Width.low),
              fill = ribbon.fill, alpha = ribbon.alpha) +
  geom_line(aes(x = Year, y = Width.mean), size = 1) +
  geom_line(aes(x = Year, y = Width)) +
  geom_point(aes(x = Year, y = Width), size = 2) +
  xlab("year") +
  ylab("growth width (mm)") +
  ylim(0, 10) +    
  theme_classic(base_family = base.family)
print(p3)

Fig. 'NenRin' growth width (by {dlm})

Un examen plus attentif montre que la taille de la zone grise, qui indique l'intervalle de confiance à 95%, est légèrement différente, mais elle est presque la même que le résultat de JAGS montré ci-dessus.

Étant donné que j'utilise habituellement Python, j'aimerais étudier l'implémentation du système de filtre Kalman de Python. J'ai essayé d'utiliser le Jupyter Notebook cette fois-ci, et j'ai réalisé qu'il donne une bonne vue de l'analyse des données de "R". Au fait, "Kalman and Bayesian Filters in Python" est publié sur Github comme exemple de Jupyter Notebook, il est donc recommandé pour ceux qui sont intéressés. (Comme il s'agit d'une assez grande quantité de contenu, je le lis petit à petit. Des animations, etc. sont également utilisées, et comme il s'agit d'un Jupyter Notebook, vous pouvez également déplacer le code directement. Le lien est décrit ci-dessous.)

Références (site Web)

• Page d'assistance "Iwanami Data Science" https://sites.google.com/site/iwanamidatascience/
• IRkernel - GitHub https://github.com/IRkernel/IRkernel --Modèle de niveau local --Logic of Blue http://logics-of-blue.com/%E3%83%AD%E3%83%BC%E3%82%AB%E3%83%AB%E3%83%AC%E3%83%99%E3%83%AB%E3%83%A2%E3%83%87%E3%83%AB/

--Reproduire "Introduction à l'analyse des séries spatio-temporelles d'état" avec R http://elsur.jpn.org/ck/

• {dlm} Manuel attaché au package (dlm: un package R pour l'analyse bayésienne des modèles linéaires dynamiques)
• "Kalman and Bayesian Filters in Python” https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python