[PYTHON] Probieren Sie das Zustandsraummodell aus (Jupyter Notebook + IR-Kernel).

Umriss des Zustandsraummodells

Das Buch ist Iwanami Data Science Vol.1. ([Besonderheit] Bayesianisches Denken und freie MCMC-Software) Das konzeptionelle Diagramm des Zustandsraummodells wird hier zitiert.

Das Zustandsraummodell hat eine Zeitreihe ** "Zustand" **, und basierend darauf wird ein ** "beobachteter Wert" ** erzeugt, der eine bestimmte Variationsmenge enthält. Das ist. Iwanami Data Science Vol.1 ist eine Besonderheit von "Bayes Inference and MCMC Free Software", aber Kapitel 3 ist ein Artikel über Zustandsraummodelle (geschrieben von Herrn Ito), und andere Kapitel erklären auch Raummodelle. Kommen Sie. In Kapitel 3 wurde erklärt, wie dies mit R und JAGS (Just Another Gibbs Sampler) gelöst werden kann.

Als ich das versuchte, versuchte ich es zuerst mit Python + PyMC3, aber es funktionierte nicht.

import pymc3
(Weggelassen)
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):   #Die folgenden zwei Zeilen sind FEHLER.Ich wollte die Zustandsgröße α definieren...
        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)
(Weggelassen)

Ich habe diese Art von Code ausprobiert, aber in der for-Schleife wird an der Stelle ein Fehler angezeigt, an der das Zustandsmengenmodell $ \ alpha_t = \ alpha_ {t-1} + \ eta $ definiert ist. Wenn Sie darüber nachdenken, ist ** pymc3.Normal () ** (** pm.Normal () ** in der Liste) ein Klassenkonstruktor zum Generieren eines Modells der Normalverteilung, also in seinem Argument Es wird als schwierig angesehen, die zu generierende Klasse (Alpha [i-1]) einzuschließen. (Ich habe das Internet unterschiedlich durchsucht, aber keine Lösung gefunden.)

Ich habe überlegt, ** PyStan ** und die Implementierung des Kalman-Filtersystems ** PyKalman ** als andere Bibliotheken zu verwenden, die solche Zustandsraummodelle (oder dynamische lineare Modelle) verarbeiten können, aber diesmal bin ich dem Buch gefolgt Ich habe mich für "R" entschieden. Das verwendete Programm ist wie folgt.

• R ver.3.2.2 --Jupyter Notebook (Notebook ist Version 4.0.6) + IRkernel ... Dieses Mal habe ich IRkernel zum ersten Mal ausprobiert. --JAGS ver.4.0.0 + {rjags} Paket
• {dlm} Paket

Informationen zur Umgebung von Jupyter Notebook + IRkernel

Jupyter Notebook ist eine aktualisierte Version des IPython Notebook, das Sie kennen. Durch Umschalten des Sprachverarbeitungssystems und des Kernels ist es möglich, andere als Python zu verarbeiten. (Es wurde nach Jupyter = Julia + Python + R benannt.) Dieses Mal habe ich IRkernel für die R-Sprache installiert und verwendet.

Das Installationsverfahren ist wie folgt.

-Installieren Sie die R-Sprache (falls erforderlich). --Installieren Sie IRkernel (https://github.com/IRkernel/IRkernel) in einer R-Umgebung.

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

--Starten Sie mit "> Jupyte Notebook" auf der Konsole. (Zuerst ist ein Fehler aufgetreten, weil einige R-abhängige Pakete fehlen. Er kann jedoch gestartet werden, indem das R-Paket gemäß der Fehlermeldung installiert wird.)

** Abb. Wenn Jupyter Notebook gestartet wird **

Wie oben gezeigt, können Sie "R" aus dem neuen Menü zur Notenerstellung auswählen.

Implementierung des Zustandsraummodells durch JAGS

Folgen Sie dem Kapitel "Modellierung von Zeitreihen und räumlichen Daten" im Buch "Iwanami Data Science Vol.1", um die folgenden Hinweise zu erstellen.

Wie bei früheren IPython-Notizbüchern können Sie Markdown-Text mit ausführbarem Code mischen. Außerdem bin ich froh, dass ich mit diesem Update von Notebook ** MathJax ** -Formeln offline rendern kann.

Ein Tipp beim Zeichnen einer Figur mit R. (Notizbuch in [1]: Dies ist der beschriebene Befehl.)

library(repr)

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

In der zweiten Hälfte dieses Hinweises zeichnen wir mit {ggplot2}. Wenn jedoch nichts angegeben ist, wird die Abbildung so generiert, dass sie der gesamten Breite des Webbrowsers entspricht. Die Zahl ist ziemlich groß (obwohl sich die meisten Menschen in einer ähnlichen Situation befinden würden), und die Balance ist nicht sehr gut, wenn sie zusammen mit Sätzen und Codes angezeigt wird. {repr} ist ein Unterpaket, das zusammen mit IRkernerl installiert wird. Mit dieser Funktion können Sie jedoch die oben genannten Optionen angeben, um die Figur in einer beliebigen Größe zu erstellen.

In diesem Hinweis wird auch der BUGS-Code beschrieben, obwohl dies nicht der Code ist, der direkt ausgeführt werden soll. (JAGS + {rjags} ist eine Spezifikation, die den BUGS-Code in einer separaten Datei vorbereitet und die Berechnung ausführt.)

Ich fand es ein großer Vorteil, Jupyter Notebook zu verwenden, um den R-Code, der das Ganze zusammenfasst, den BUGS-Code, der das Simulationsmodell definiert, die Ausgabe der Berechnungsergebnisse und die Informationen, die dazu neigen, getrennt zu streuen, kombinieren zu können.

Der auf diese Weise vom lokalen Modell berechnete Übergang der jährlichen Ringbreite von Sugi in Kyoto wird wie folgt dargestellt.

Fig. 'NenRin' growth width (by JAGS)

Die gestrichelte Linie mit dem Marker ist der beobachtete Wert, die dicke durchgezogene Linie ist der durch die Berechnung von MCMC geschätzte Wert und die Grauzone ist der Kreditabschnitt von 95%. (Der R-Code ist auf der Support-Website von "Iwanami Data Science" erhältlich.)

Bestätigung durch {dlm} Paket

JAGS ist ein Implementierungssystem von MCMC (Marcov-Kette Monte Carlo). Zusätzlich zu MCMC gibt es eine Implementierung des Kalman-Filtersystems als Berechnungsmethode, die Zustandsraummodelle (dynamische lineare Modelle) verarbeiten kann. In "R" werden jedoch häufig die Pakete {dlm} und {KFAS} verwendet. Diesmal habe ich das {dlm} -Paket ausprobiert.

Der Unterschied zwischen MCMC- und Kalman-Filter wird in "Iwanami Data Science" ausführlich erläutert. (Ich werde zitieren ...) "Der Unterschied zwischen dem Kalman-Filter und MCMC besteht darin, dass die periphere posteriore Verteilung für den Zustand $ x_i $ berechnet werden kann, aber es sieht aus wie $ s ^ 2 $, $ \ sigma ^ 2 $. Für die Parameter des Modells wird die Schätzung (Punktschätzung) durchgeführt, wobei die Streuung der posterioren Verteilung ignoriert wird. “Die Streuung der Parameter der Dispersionsrelation, sozusagen die„ Variation der Variation “, ist (persönlich) nicht sehr interessiert. ， Dieser Punkt spielt zunächst keine Rolle.

Unten wird der R-Code aus dem erstellten Notizbuch zitiert.

## 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)

Mit dem Obigen kann das Berechnungsergebnis des Kalman-Filters erhalten werden. Zeichnen Sie dieses Ergebnis mit {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})

Obwohl die Größe des grauen Bereichs, der das 95% -Konfidenzintervall anzeigt, im Detail leicht unterschiedlich ist, entspricht sie fast dem oben gezeigten Ergebnis von JAGS.

Da ich normalerweise Python verwende, möchte ich die Implementierung des Kalman-Filtersystems von Python untersuchen. Ich habe diesmal versucht, das Jupyter-Notizbuch zu verwenden, und festgestellt, dass es einen guten Überblick über die Datenanalyse von "R" gibt. Übrigens wird "Kalman and Bayesian Filters in Python" auf Github als Beispiel für Jupyter Notebook veröffentlicht, daher wird es allen Interessierten empfohlen. (Da es sich um eine ziemlich große Menge an Inhalten handelt, lese ich sie nach und nach. Es werden auch Animationen usw. verwendet. Da es sich um ein Jupyter-Notizbuch handelt, können Sie den Code direkt verschieben. Der Link wird unten beschrieben.)

Referenzen (Website)

• Support-Seite "Iwanami Data Science" https://sites.google.com/site/iwanamidatascience/

• IRkernel - GitHub https://github.com/IRkernel/IRkernel

• Modell auf lokaler Ebene - Logik von Blau 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/

• Wiedereinführung von "Einführung in die Analyse von Zustandsraum-Zeitreihen" mit R http://elsur.jpn.org/ck/

• {dlm} Dem Paket beigefügtes Handbuch (dlm: ein R-Paket für die Bayes'sche Analyse dynamischer linearer Modelle)

• "Kalman and Bayesian Filters in Python” https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python