ブログ

ορεσικα ψομαναι κυσο βλογ

Rによる時系列モデリング入門 第1章 データのグラフ化まで

この文書は

2020年12月11日に岩波書店から出版された、北川源四郎著の『Rによる時系列モデリング入門』について、最初のデータの読み込みから基本的なグラフ化までを確認のために出力してみたものである。やってみて、いくつか気づいた点があるのでメモを残す。

必要なパッケージをロード

本書で必要な関数とデータを全て含むTSSSパッケージをロードする。

library(TSSS)

以下のパッケージはRを用いた現代的なデータ解析には必要不可欠である。

library(tidyverse)

以下のパッケージはtsクラスの変数をデータフレームに変換するために用いる。

library(ggpmisc)

データ取り込み

data("HAKUSAN")
data("Sunspot")
data("Temperature")
data("BLSALLFOOD")
data("WHARD")
data("MYE1F")

以下の3つのデータセットはTSSSパッケージに含まれていないようである。

  • Nikkei225
  • rainfall
  • Haibara

webで検索したら著者が公開しているwebページが見つかった。 恐らく以下のこれらが最新版の元データであろう。

Nikkei225 <- read_csv("http://www.mi.u-tokyo.ac.jp/mds-oudan/lecture_document_2019_math7/%E6%99%82%E7%B3%BB%E5%88%97%E3%83%87%E3%83%BC%E3%82%BF/nikkei225_new.csv")

Rainfall <- read_csv("http://www.mi.u-tokyo.ac.jp/mds-oudan/lecture_document_2019_math7/%E6%99%82%E7%B3%BB%E5%88%97%E3%83%87%E3%83%BC%E3%82%BF/rainfall_new.csv")

Haibara <- read_csv("http://www.mi.u-tokyo.ac.jp/mds-oudan/lecture_document_2019_math7/%E6%99%82%E7%B3%BB%E5%88%97%E3%83%87%E3%83%BC%E3%82%BF/haibara_new.csv")

ただし、これらのURLがいつまで有効かは不明である。

CSVデータの取り込み例は省略する。上のread_csvによるwebからのデータ読み込みで例示は十分だろう。

plot例

本書では素朴なplot関数によるグラフ化が説明されているが、より実践的なggplotを用いてグラフ化を行う。

ggplot(data = HAKUSAN) +
  aes(x = as.numeric(row.names(HAKUSAN)), y = YawRate) +
  xlab("(a) 船舶の方向角速度") + ylab("") +
  geom_line()

f:id:XX-Prime:20201215002828p:plain

Sunspot.df <- try_data_frame(Sunspot)
ggplot(data = Sunspot.df) +
  aes(x = time, y = Sunspot) +
  xlab("(b) 太陽黒点数") + ylab("") +
  geom_line()

f:id:XX-Prime:20201215002831p:plain

Temperature.df <- data.frame(number = Temperature)
ggplot(data = Temperature.df) +
  aes(x = as.numeric(row.names(Temperature.df)), y = number) +
  xlab("(c) 東京の日最高気温") + ylab("") +
  geom_line()

f:id:XX-Prime:20201215002834p:plain

BLSALLFOOD.df <- try_data_frame(BLSALLFOOD)
ggplot(data = BLSALLFOOD.df) +
  aes(x = time, y = BLSALLFOOD) +
  xlab("(d) 食品産業従事者数") + ylab("") +
  geom_line()

f:id:XX-Prime:20201215002837p:plain

WHARD.df <- try_data_frame(WHARD)
ggplot(data = WHARD.df) +
  aes(x = time, y = WHARD) +
  xlab("(e) 卸売高") + ylab("") +
  geom_line()

f:id:XX-Prime:20201215002841p:plain

MYE1F.df <- data.frame(number = MYE1F)
ggplot(data = MYE1F.df) +
  aes(x = as.numeric(row.names(MYE1F.df)), y = number) +
  xlab("(f) 地震波(東西方向)") + ylab("") +
  geom_line()

f:id:XX-Prime:20201215002845p:plain

ggplot(data = Nikkei225) +
  aes(x = as.numeric(row.names(Nikkei225)), y = nikkei225) +
  xlab("(g) 日経 225 平均株価") + ylab("") +
  geom_line()

f:id:XX-Prime:20201215002849p:plain

Haibara.tidy <- Haibara %>% 
  gather(key = measure, value = value, 地下水位, 気圧)
ggplot(Haibara.tidy, aes(x = as.numeric(日時), y=value)) + 
  geom_line() + facet_grid(measure ~ ., scales = "free") + 
  xlab("(h) 地下水位(下段)と気圧(上段)") + ylab("")

f:id:XX-Prime:20201215002852p:plain

HAKUSAN.tidy <- HAKUSAN %>% 
  mutate(index = as.numeric(row.names(.))) %>% 
  gather(key = measure, value = value, Rolling, Pitching, Rudder)
ggplot(HAKUSAN.tidy, aes(x = index, y=value)) + 
  geom_line() + facet_grid(measure ~ ., scales = "free") + 
  xlab("(i) 船舶の縦揺れ(上段), 横揺れ(中段), 舵角(下段)") + ylab("")

f:id:XX-Prime:20201215002855p:plain

気づきなど

tsクラスのデータはそのままではggplotでグラフ化できない

(b)(d)(e)をグラフ化する際に、ggpmiscのtry_data_frame()でデータフレームに変換してからggplotでグラフ化した。

参考

https://www.karada-good.net/analyticsr/r-453

facetで順序を制御するのはどうやら面倒

(h)(i)をグラフ化する際に、もとのデータフレームをtidyな形式に変換した後にfacet_gridで縦に並べて表示したが、順番はアルファベット順になってしまうようで本の順番通りには配置できなかった。ファクター変数でやや無理やり順番を制御することもできるようだが、面倒なのでggplotの表示順をそのまま採用した。

(h)地下水位のデータは長い欠損期間があるようだ

本では行数をx軸のインデックスとして使っているようだが、地下水位データには「日時」という列があり、今回はこれをx軸として図示した所、図のように長い欠損期間が二ヶ所あることが分かった。

本文ではこの長い欠損期間については特に触れられていないようである。

この文書をどうやって作成したか

rmarkdownで作成したが、出力された文書をHatena Blogに載せる方法がわからなかった。以下の解説を参考にmdファイルを生成してコピペした。画像の添付が面倒である。どうにかならんもんか。

参考

mfavoritey.hatenablog.com

M2X無料アカウントが停止するのでデータをエクスポートする

M2X無料アカウント停止

メールを見逃していたが、日本時間で2020/04/16 02:00 にM2Xの無料アカウントが停止される。2020/02/15にそのメールが着ていたが、件名からはそのような重大なメールであることがわからず読んでいなかった。

データエクスポート

  1. データエクスポートコマンド発行→ジョブID発行
  2. ジョブIDでジョブ状態を確認
  3. ジョブ終了でデータエクスポート用URLが表示されるのでwebブラウザなどでダウンロード

という流れになる。

適当なLinuxホストにログインして(今回はWSLを利用)データエクスポートコマンドをcurlで実行する。devices/以降はデバイスID、X-M2X-KEY:にはデバイスAPIキーを指定する。

curl -i "https://api-m2x.att.com/v2/devices/<DEVICE-ID>/values/export.csv" -H "X-M2X-KEY: <DEVICE-API-KEY>"

次のような応答でエクスポートジョブが受理されたことがわかる。

HTTP/1.1 202 Accepted
Content-Type: application/json; charset=UTF-8
Location: https://api-m2x.att.com/v2/jobs/<JOB-ID>
X-M2X-VERSION: v2.112.2
Vary: Accept-Encoding
X-RateLimit-Limit: 25
X-RateLimit-Remaining: 24
X-RateLimit-Reset: 1586865919
Content-Length: 21

{"status":"accepted"}

最後のacceptedで受理がわかり、ジョブ確認用のURLがLocation:の次に表示されている。

エクスポートジョブの状態を以下のコマンドで確認する。ここで指定するjobs/以降のIDは上のジョブ確認用のIDで、X-M2X-KEY:にはデバイスキーではなくマスターAPIキーを指定する。マスターAPIキーはアカウント設定からMaster Keyで確認できる。

curl -i "https://api-m2x.att.com/v2/jobs/<JOB-ID>" -H "X-M2X-KEY: <MASTER-API-KEY>"
HTTP/1.1 200 OK
Content-Type: application/json; charset=UTF-8
X-M2X-VERSION: v2.112.2
Vary: Accept-Encoding
X-RateLimit-Limit: 25
X-RateLimit-Remaining: 24
X-RateLimit-Reset: 1586866332
Content-Length: 220

{"id":"<job-id>","state":"working","output":null,"errors":null,"started":"2020-04-14T12:05:18.160Z","finished":null,"created":"2020-04-14T12:05:17.996Z","updated":"2020-04-14T12:05:18.160Z"}

まだstateはworkingなのでエクスポートジョブが実行中である。データサイズにもよるが、数分後にもう一度ジョブ確認のコマンドを発行すると、今度は以下のような応答が得られる。

{"id":"<job-id>","state":"complete","output":null,"errors":null,"started":"2020-04-14T12:05:18.160Z","finished":"2020-04-14T12:23:13.818Z","created":"2020-04-14T12:05:17.996Z","updated":"2020-04-14T12:23:13.819Z","result":{"url":"http://export.m2x.sl.attcompute.com/<DEVICE-ID>%20-%20<DEVICE-NAME>%20-%202020-04-14T12:22:53.270Z.csv"}}

url:以降に、csvフォーマットのファイルダウンロードURLが書かれているので、これをコピーしてwebブラウザなどでダウンロードする。もちろんcurlでダウンロードしても良いが、webブラウザのほうが進捗がわかりやすい。

M2Xからの移行先は

どこがいいだろうか?

買い占め行動(買いあさり)

原書房「民間防衛」
35ページ

一たび輸入がとまったとなれば、買いあさりが行われやすい。
この買いあさりは、すべての人が充分な資金力を持っているわけではないから、
反社会的な行動として、これを非難しないわけにはいかない
さらにまた、買いあさりが起これば、流通機構に平常以上の負担がかかるので、
日常生活用品の供給が困難になり、ひいては国内の在庫品の分配が不公平になる
このような買いあさりに加わることは、自分の社会連帯意識の低さ、
つまり、自分勝手なことをしたのでは社会が成り立たないという意識が欠けていることを、
暴露するだけでなく、自分が有事の際の備え、貯えを怠っていたことを証明するにほかならない

台風による気圧変化の記録

台風19号による気圧変化を二箇所で比較する

BME680で二箇所の環境データをM2Xに記録し続けている。今回の台風19号について気圧変化をメモしておく。

埼玉県鴻巣市の気圧グラフ

台風の目がかなり近くを通過したので、気圧変化は尖点を持つ特徴的なものになった。

f:id:XX-Prime:20191013173650p:plain
埼玉県鴻巣市 台風19号 気圧変化

最低気圧は971.9hPaであった。

滋賀県大津市の気圧グラフ

目はだいぶ遠くを通ったので、尖点を持たないグラフになった。

f:id:XX-Prime:20191013174108p:plain
滋賀県大津市 台風19号 気圧変化

最低気圧は976.4hPaであった。

二箇所のグラフの時間的なズレを見ると台風のだいたいの速度を概算できる。面倒なんでやらないけど。

switchbot カーテン

カーテン自動開閉マッシーン

switchbotがまた面白い製品を出すようだ。

ks.switch-bot.com

段差も越えて動く

f:id:XX-Prime:20191013205305g:plain
rail moving

日本で一般的な角型カーテンレールで動作するか?

たぶん問題ない。

カーテン以外の用途

レールなら割と何でもよさそうなので、別にカーテンだけに用途を絞る必要もなかろう。1次元的な動作を必要とする、自動化したい用途に一般的に使えるはずだ。本体にフックのような構造があると色々なモノを引っ掛けられるので便利だが、画像を見る限りそのようなものは無いので自分で何とかすればOKだろう。ある程度パワーがあるなら、2つ組み合わせて洗濯物を自動的に取り込む機構もかんたんに作れそう。

まとめ

買う。カーテン以外の用途を見つけたい。IoTのデモとかで、コレをつかっておバカで派手なモノを作る輩が現れるのが目に浮かぶ・・・。

M2Xの調子が悪い

珍しくM2Xのサーバが落ちている

f:id:XX-Prime:20191005131241p:plain
M2X_down
ほとんど落ちないM2Xのサーバが珍しく停止しているのでメモ。

M2X状態モニタを見る

上の画像でリンクがある、M2X状態モニタにアクセスしてみると、現在の状態とインシデントログを見ることができる。

f:id:XX-Prime:20191005131953p:plain
status-m2x.att.com
こういうながーいスクリーンキャプチャは、chromeのデフォルト機能でできることを今知った。参考リンク
requlog.com

2019年10月3日に障害が発生して一旦収まったものの、問題は悪化しているようだ。
現時点でレスポンスタイムが0に見えるのは、単にサーバが機能していないからである。

まとめ

珍しいこともあるもんだ。
複数サイトへのデータアップロードをしておけばよかった。
今はどのあたりのサイトがいいんだろうか。
M2Xが落ちるという想定はしていなかった。