[PYTHON] Vorhersage der Zukunft der Bevölkerungsentwicklung in Numazu City durch Zeitreihen-Regressionsanalyse mit Prophet

Was du machen willst

・ Prognostizieren Sie die zukünftige Bevölkerung von Numazu City basierend auf der Demografie der letzten Jahre mit Pythons Zeitreihen-Vorhersagebibliothek "Facebook Prophet". ・ Analysieren Sie die Analyseergebnisse. ⇒ Abstimmung ⇒ Erneut analysieren

Umgebung

**　Windows10 64bit** **　Anaconda3** Virtuelle Umgebung: Python (3.7.3) Bibliotheken: fbprophet (0.6), pandas (1.1.2), matplotlib (3.0.3)

Nutzungsdaten

Liste "Bevölkerung und Haushalt" der Stadt Numazu

• Verwenden Sie die aggregierte monatliche Grundgesamtheit aus den jährlichen Dateien auf der oben genannten Site als Eingabedatei ** [dataset.csv] **. ・ Ursprünglich enthält die offenen Daten von Numazu City einen demografischen Bericht. Offener Datenkatalog von Fuji no Kuni (Stadt Numazu) Dies sind jährliche Statistiken, aber da wir monatlich als Granularität der Daten für die Zeitreihenanalyse angenommen haben, Dieses Mal werden wir die Informationen im Basic Resident Register auf der offiziellen Website der Stadt, die monatliche Daten enthält, bequem verwenden. -Der Preis ist, dass Sie das PDF in eine Tabelle formatieren müssen. Es war mühsam. ・ Es besteht aus zwei Reihen von Jahr-, Monats- und Haushaltspersonal. Da dies ein Versuch ist, gibt es keine externen Variablen. ・ 2013 / 04-2020 / 08 werden die Lerndaten sein.

dataset.csv

Vorausschauende Analyse

Bibliothek importieren

import pandas as pd
import matplotlib.pyplot as plt
from fbprophet import Prophet

Datenerfassung

Jinkoudata = pd.read_csv("dataset.csv")
Jinkoudata.columns = ['ds','y']

-Importieren Sie Bevölkerungsdaten in den Datenrahmen [Jinkou-Daten].

• Ändern Sie die Spaltennamen [Jahr / Monat] und [Haushaltspersonal (insgesamt)] in [ds] [y](gemäß der Namenskonvention des Propheten).

Modellgenerierung

model = Prophet(weekly_seasonality=False,daily_seasonality=False)

・ Da die Analyse dieses Mal monatlich durchgeführt wird, werden Parameter angegeben, die den Wochenzyklus und den Tageszyklus nicht berücksichtigen.

• Außerdem können die folgenden Parameter beim Erstellen einer Instanz beliebig angegeben werden. Parameterwerte in []. Der als Beispiel angegebene Wert ist der Standardwert.

●growth='linear' [logistic/linear] Geben Sie das Analysemodell als linear oder logistisch an.

●changepoints=None [DATE(YYYY-MM-DD)] Beliebige Bezeichnung von Änderungspunkten von Trends (langfristige Trendschwankungen). Im Propheten werden Änderungspunkte automatisch extrahiert und in der Vorhersage wiedergegeben, aber anstatt die automatische Erkennung von Änderungspunkten nach Bedarf zu verwenden, Es wird mit der Absicht verwendet, die Position potenzieller Änderungspunkte mithilfe dieses Parameters manuell anzugeben.

** ● n_changepoints ** = 25 [beliebige Konstante] Geben Sie die Anzahl der zu erkennenden Änderungspunktkandidaten an. Je größer der Wert, desto einfacher ist es, viele Änderungspunkte zu erkennen. Wenn dies zu groß gemacht wird, werden Änderungspunkte möglicherweise übermäßig erkannt und die Vorhersagegenauigkeit kann abnehmen.

** ● jährliche_Saisonalität ** = 'auto' [true / false / auto / beliebige Konstante] Geben Sie die Periodizität (Jahr) an. Der Grad der Berücksichtigung kann durch Angabe einer beliebigen Zahl angepasst werden (10 ist der angegebene Wert, wenn nicht angegeben).

** ● wöchentliche Saisonalität ** = 'auto' [true / false / auto / beliebige Konstante] Geben Sie die Periodizität (Woche) an. Der Grad der Berücksichtigung kann durch Angabe einer beliebigen Zahl angepasst werden (10 ist der angegebene Wert, wenn nicht angegeben).

** ● daily_seasonality ** = 'auto' [true / false / auto / beliebige Konstante] Geben Sie die Periodizität (Tag) an. Der Grad der Berücksichtigung kann durch Angabe einer beliebigen Zahl angepasst werden (10 ist der angegebene Wert, wenn nicht angegeben).

●holidays=None [DataFrame] Geben Sie den Datenrahmen für das Ereignisdatum an. Integrieren Sie Feiertage und Daten mit Ereignissen in die Analyse und spiegeln Sie diese in der Prognose wider.

** ● changepoint_prior_scale ** = 0,05 [beliebige Konstante] Geben Sie die Streuung der Laplace-Verteilung an. Je größer dieser Wert ist, desto einfacher ist es, den Änderungspunkt zu erkennen. Es kann jedoch zu flexibel sein. Gleich wie n_changepoints.

** ● changepoint_range ** = 0,8 [beliebige Konstante] Geben Sie das Datenverhältnis zum Erkennen von Änderungspunkten an. Da der angegebene Wert 0,8 ist, werden die Änderungspunkte zu 80% aus dem Anfang der Trainingsdaten extrahiert. Wenn sich der Trend am Ende der Daten beispielsweise stark bewegt, setzen Sie ihn auf "1.0". Der Änderungspunkt wird anhand des gesamten Datenbereichs einschließlich der Schwankung in der Endstufe erfasst.

Es scheint verschiedene andere Dinge zu geben, aber ich werde es weglassen, weil das Obige das einzige ist, was mit der Menge an Wissen bisher berücksichtigt werden kann.

Lernen

model.fit(Jinkoudata)

fit() ・ Lernen Sie, indem Sie die Trainingsdaten mit Fit angeben.

Datenrahmengenerierung

future = model.make_future_dataframe(periods=16,freq='M')

make_future_dataframe()

• Erstellen Sie eine zukünftige Zeitreihenzeile, indem Sie sie dem Datenrahmen hinzufügen, der als Trainingsdaten für das erstellte Modell verwendet wird. ・ Hier zusätzlich zum angepassten Datenrahmen (Jinkoudata: 2013 / 04-2020 / 08), Wir haben Datenrahmen für die nächsten 16 Monate (2020 / 09-2021 / 12) erstellt. ・ Infolgedessen werden [Zukunft] Daten mit einem Zeitraum von 2013/04 bis 2021/12 sein.

future

Prognose

predict = model.predict(future)

• Geben Sie den zuvor erstellten Zeitreihendatenrahmen als Argument an und sagen Sie ihn voraus. -Die vorherzusagende Zeitreihe ist 2013/04 bis 2021/12. Dies ist die Zeitreihe, die in Zukunft gespeichert wird.

** Vorhersageergebnis [vorhersagen] **

Es gibt andere Informationen, die vorhergesagt und im Datenrahmen gespeichert werden. Da es schwierig ist, den numerischen Wert zu erfassen, werden wir hier die Elemente überprüfen, die sich auf den vorhergesagten Wert beziehen. ・ Yhat: Voraussichtlicher Wert ・ Yhat_lower: Voraussichtliche Untergrenze ・ Yhat_upper: Voraussichtliche Obergrenze "Prognostizierter Wert", "Prognostizierter unterer Grenzwert" und "Prognostizierter oberer Grenzwert" werden in den obigen Spalten gespeichert. "Vorhersagewert" ist ein sehr wahrscheinliches Vorhersageergebnis, und "vorhergesagte untere Grenze" und "vorhergesagte obere Grenze" sind die unteren und oberen Grenzen des Schwankungsbereichs des vorhergesagten Werts.

Analyse

model.plot(predict)

• Zeichnen Sie die Vorhersage mit [plot]. ・ Schwarze Punkte sind Trainingsdaten, blaue Linien sind vorhergesagte Werte (yhat) und blaue Bänder sind vorhergesagte Schwankungsbreiten (yhat_lower: yhat_upper).

** [Was können Sie hier lesen] ** ■ Aufgrund der Prognose ist der Trend des Rückgangs im Jahr 2021 im gleichen Tempo wie zuvor (etwa 1000 Personen oder weniger / Jahr). ■ Bis Ende 2021 das Niveau von 193.000 durchbrechen ■ Die Genauigkeit der Vorhersage ist hoch, da der vorhergesagte Wert und das tatsächliche Ergebnis nicht getrennt sind (vielmehr sind die Daten einfach). ■ Es ist möglich, Vorhersagen mit einem relativ kleinen Schwankungsbereich der vorhergesagten Werte bis etwa ein halbes Jahr im Voraus zu treffen. (Die Lerndaten haben von 2021/03 bis 2020/08 rapide zugenommen.)

model.plot_components(predict)

-Das Zeichnen für jeden Faktor wie Trendschwankungen und saisonale Schwankungen erfolgt mit [plot_components]. ・ Der [Trend] -Term ist der Trend (langfristige Trendschwankung). Ich zeichne, welche Art von Fluktuation es als jährlicher Übergang ist. ・ Die [jährliche] Laufzeit ändert sich saisonal. Es zeigt die Schwankungen während des Jahres. ・ Trends zeigen in kurzer Zeit keine starke Zunahme oder Abnahme und beobachten eine jährliche Zunahme oder Abnahme der Skala, die saisonale Schwankungen wiederholt.

** [Was können Sie hier lesen] ** ■ In den letzten 8 Jahren hat es eine schöne Linienform angenommen. Nach dem Jahr des Lernens (2013) gibt es keine Anzeichen für eine jährliche Erholung ■ Als Schwankung innerhalb des Jahres nimmt sie im Wesentlichen weiter ab, wird jedoch im Oktober und November leicht zunehmen. (Da saisonale Schwankungen analysiert und durch Wellenformen ausgedrückt werden, werden monatliche Daten wie diese Zeit verwendet Beachten Sie, dass erwartet wird, dass die Wellenform größer wird, wenn das Datendatumsintervall größer ist.) (Obwohl die tatsächlichen Daten kaum zunehmen, ist die Abnahmerate im Oktober und November im Vergleich zu anderen Daten erheblich zurückgegangen. Da es in den letzten Jahren einen leichten Anstieg gezeigt hat, scheint es einen Ort zu geben, an dem es sich als Vorhersage zu einem Anstieg entwickelt hat.) ■ Der Abwärtstrend ist leicht allmählich (Obwohl der Rückgang jedes Jahr im März stark ist, ist auch die Abnahmerate selbst rückläufig. Zum Beispiel ist der Rückgang im März selbst auf Beschäftigungs- / Aufstiegsereignisse zurückzuführen. Es kann spekuliert werden, dass die Verringerung der Abnahmerate auf die sinkende Geburtenrate zurückzuführen ist. Um die Details zu erfahren, soziale Dynamik wie Ein- und Ausziehen, wie die Anzahl der Geburten und Todesfälle Offene Daten wie natürliche Dynamik sind erforderlich.

Auswertung

Es muss bewertet werden, wie viel Vorhersagegenauigkeit das erstellte Modell aufweist. Dieses Mal werden wir mithilfe der Kreuzungsüberprüfung bewerten.

Gegenüberstellung

from fbprophet.diagnostics import cross_validation
df_cv = cross_validation(model, initial='1095 days', period='182 days', horizon = '365 days')

cross_validation() ・ Anfänglich: Wann soll mit der Prognose begonnen werden? ・ Zeitraum: Zeitraum zum Verschieben von Daten ・ Horizontal: Erwarteter Zeitraum Im Fall des obigen Deklarationsbeispiels ab dem 1095. Tag (Initiale) 1 bis 1094 Tage Daten für den Trainingsdatenzeitraum (3-facher Horizont, sofern nicht anders angegeben), Prognostizieren Sie den Bewertungszeitraum vom 1095. Tag bis zu 365 Tagen (Horizont). Verschieben Sie nach der Analyse der Prognose den Startpunkt des Prognosezeitraums um 182 Tage (Zeitraum) zurück. Wir prognostizieren 365 Tage ab dem 1277. Tag als nächsten Bewertungszeitraum. Dies wird bis zum Ende der Daten wiederholt und das Vorhersageergebnis gespeichert. Dieses Vorhersageergebnis wird als Bewertungsergebnis verwendet, indem es durch eine Skala ersetzt wird, die durch das später beschriebene Verfahren leicht bestätigt werden kann.

** Kreuzverifizierungsergebnis [df_cv] **

Bewertung der Bewertungsskala

from fbprophet.diagnostics import performance_metrics
df_p = performance_metrics(df_cv)

performance_metrics() Berechnet die Bewertungsskala für das Ergebnis der Kreuzungsüberprüfung [df_cv] und speichert das Ergebnis.

** Berechnungsergebnis der Bewertungsskala [df_p] **

Die Bewertungsskala ist eine Berechnung der Ergebnisse für jeden Bewertungsindex, und die Typen sind wie folgt.

** ● MSE [mittlerer quadratischer Fehler] ** Der Durchschnittswert des Quadrats des Absolutwerts, der die Differenz zwischen dem tatsächlichen Wert und dem vorhergesagten Wert darstellt.

** ● RMSE [quadratischer Mittelwertfehler] ** Die Quadratwurzel des Quadrats der Differenz zwischen dem tatsächlichen Wert und dem vorhergesagten Wert.

** ● MAE [mittlerer absoluter Fehler] ** Der Durchschnitt der absoluten Werte der Differenz zwischen dem tatsächlichen Wert und dem vorhergesagten Wert.

** ● MAPE [mittlerer absoluter prozentualer Fehler] ** Durchschnitt des Verhältnisses der Differenz zwischen dem tatsächlichen Wert und dem vorhergesagten Wert zum tatsächlichen Wert.

MSE, RMSE und MAE werden unterschiedlich berechnet, aber die Abweichungen in der Vorhersage werden im Durchschnitt basierend auf tatsächlichen Messungen (diesmal Population) ausgedrückt. MAPE drückt die Abweichung der Vorhersage in Prozent (%) aus. Da die diesmal behandelten Zahlen groß sind, möchte ich die Bewertungsergebnisse mit MAPE (Average Absolute Percentage Error) bestätigen, das intuitiv leicht zu verstehen ist.

Visualisierung der Bewertungsskala durch MAPE (durchschnittlicher absoluter prozentualer Fehler)

from fbprophet.plot import plot_cross_validation_metric
fig = plot_cross_validation_metric(df_cv, metric='mape')

plot_cross_validation_metric() Sie können die im vorherigen Abschnitt berechnete Bewertungsskala visualisieren. Mit dem in der Kreuzungsüberprüfung als X-Achse angegebenen Horizont (Auswertungszeitraum) Der MAPE-Wert für die Anzahl der ausgeführten Auswertungen wird als Punkt gezeichnet. Der Durchschnitt wird als blaue Linie angezeigt.

** [Was können Sie hier lesen] ** ■ Die Genauigkeit der Gesamtvorhersage ist nicht gering ⇒ Selbst wenn Sie den maximalen Abweichungswert betrachten, beträgt der Mape-Wert etwa 0,4%.

Wenn Sie dem Durchschnitt folgen, ändert sich der Fehler zwischen 0,05 und 0,15%, was ungefähr 100 bis 300 Personen entspricht. Da es sich um eine Vorhersage bis zum letzten handelt, gibt es keinen Index, der korrekt sein könnte. Da es jedoch viele Zeiträume gibt, in denen sich die tatsächlichen Messdaten monatlich um 100 Personen bewegen, besteht der Unterschied darin, dass sie in einem Aggregationszeitraum ausgefüllt werden können, sodass es sich nicht um einen großen Fehler handelt. Das können Sie beurteilen. ■ Die Genauigkeit nimmt mit Ablauf der Periode ab dem Beginn der Vorhersage ab ⇒ Natürlich nimmt die Genauigkeit im Laufe der Monate tendenziell ab (die Verteilung der Mape nimmt zu).

Tuning

Raum zur Verbesserung der Modellgenauigkeit

Da das erstellte Modell nicht immer das beste ist, kann die Genauigkeit möglicherweise verbessert werden, indem Vorhersagen unter Verwendung des abgestimmten Modells getroffen werden. Hier werde ich auch den Prozess der Beurteilung schreiben: "Gibt es Raum zum Stimmen?"

Abstimmungsperspektive

Bei der Entscheidung, ob eine Abstimmung durchgeführt werden soll oder nicht, muss die Art der für das Training verwendeten Daten berücksichtigt werden. Dieses Mal werden wir aus den folgenden Perspektiven betrachten.

** ● Gesetzliche Saisonalität / Reitsaisonalität ** Schwankungen der tatsächlichen Daten sind entweder additiv (der Bereich und die Größe der Schwankungen sind unabhängig von der Größe der Zielvariablen konstant) oder multipliziert (je größer die Zielvariable, desto größer der Schwankungsbereich). Oder

** ● Anzahl der extrahierten Änderungspunkte ** Ist die Anzahl der Änderungspunkte in den tatsächlichen Daten angemessen?

** ● Extraktionsbereich der Änderungspunkte ** Ist der Extraktionsbereich des Änderungspunktes der tatsächlichen Daten angemessen?

Beurteilung der Abstimmung auf dieses Modell

** ● Gesetzliche Saisonalität / Reitsaisonalität ** Die Zielvariable sind diesmal Bevölkerungsdaten. "Tod" und "Geburt", die als einer der Faktoren für die Zunahme und Abnahme der Bevölkerung angesehen werden, sollen alt und jung sein Weil angenommen wird, dass sich das Ausmaß der Fluktuation je nach Bevölkerungszahl ändern wird. ** Wir beurteilen, dass diese Daten einen leichten Multiplikationsaspekt haben **. (Obwohl die Gesamtzahl von Jahr zu Jahr von der Darstellung des Vorhersageergebnisses abnimmt, wird sie von dieser begleitet. Sie sehen, dass der jährliche Rückgang gering ist, wenn auch nur geringfügig. (Abgesehen davon werden "Einzug" und "Auszug" auch von der Bevölkerungsbevölkerung beeinflusst. In diesen Fällen sind auch die Anzahl und der Umfang der Unternehmen / Schulen, die Empfänger des Ein- und Auszugs sind, stark betroffen. Es war schwierig, nur mit Bevölkerungsdaten abzuschließen. Die Modellierung von Propheten kann mit dem Parameter [Seasonality_Mode] oder "Additive Seasonality" erfolgen. Sie können "legale Saisonalität" angeben. Wenn dies jedoch weggelassen wird, wird das Modell standardmäßig mit "additive Saisonalität" erstellt. Es gibt Raum, um zu versuchen, ein Modell mit "legaler Saisonalität" als Abstimmung zu erstellen. ⇒ Ziel </ font> optimieren

#Modellgenerierung (legale Saisonalität)
#Generierungszeitparameter hinzugefügt[seasonality_mode = "multiplicative"]
model_multi = Prophet(weekly_seasonality=False,daily_seasonality=False,seasonality_mode = "multiplicative")
#Lernen
model_multi.fit(Jinkoudata)
#Prognose
predict_multi = model_multi.predict(future)
#Vorhersageergebniszeichnung
fig = model_multi.plot(predict_multi)

#Gegenüberstellung
df_cv_multi = cross_validation(model_multi, initial='1095 days', period='182 days', horizon = '365 days')
#Bewertung der Bewertungsskala
df_p_multi = performance_metrics(df_cv_multi)
#Visualisierung der Bewertungsskala durch MAPE (durchschnittlicher absoluter prozentualer Fehler)
fig_multi = plot_cross_validation_metric(df_cv_multi, metric='mape')

#Einstellung der Tabellengröße
plt.figure(figsize=(12,6))
#Einstellung der Achsenrichtung
plt.xlabel('Year month')
plt.ylabel('Jinkou')
#Zeichnungsbereich angeben
plt.xlim(82, 105)
plt.ylim(190000, 196000)
#Ändern Sie die Skalenanzeige
plt.xticks([82,88,94,100,105], ["2020-01","2020-07","2021-01","2021-07","2021-12"])
#Zeichnung jedes numerischen Wertes
plt.plot(Jinkoudata.y,label = 'Jinkoudata')
plt.plot(predict.yhat,label = 'Predict_add')
plt.plot(predict_multi.yhat,label = 'Predict_multi')
plt.legend()