R言語 (2)

Time-Series

さて、前回の続きで、血圧測定データの分析。2012年からのデータが有りますからね。これを時系列的に分析して、劣化具合を調べてみましょって訳。

self data

過去にこの血圧データを分析するのを色々な言語でやったぞ。scheme,python,haskell,go-lang,awk等だ。そして今回はR言語って訳。自分の分身データだと気合が入るってもんです。

前回のスクリプトを少し整理した。

bld <- read.csv("bld.csv", header = F)
colnames(bld) <- c("yy","mm","dd","hh","hi","lo","pl")
am  <- subset(bld, bld$hh < 12)

m.am <- aggregate(cbind(hi,lo,pl) ~ mm + yy, data=am, mean)

まず、元データのbld.csv内の年の扱いが西暦下2桁だったのを4桁に改めた。sedを少し改変して年の頭に"20"を追加するだけ。ヘッダーは無し。だから、Rスクリプト中で読み込む時、header = F を指定。 FALSEが正式名称だけど、省略してます。

ヘッダーに変わる列名を colnamedsで追加してます。unix文化では、データに余計な物を付けないんで、受け取った方で何とかするのが正しい方法だと思うよ。

月毎に平均を取る方法、前回から少し変えて、最高血圧、最低血圧、脈拍も同時に計算させてます。更に、グルーピングの順番を月、年の順にした。こうしておけば、後で並べ替えは必要なくなります。

> head(m.am, n=3)
  mm   yy       hi       lo       pl
1  1 2012 125.0645 79.32258 59.29032
2  2 2012 123.2759 79.48276 59.13793
3  3 2012 120.9667 78.43333 61.83333

副作用が一つあって、列順が、月年になってしまっています。まあ欧米風と笑って許しましょう。それから、R語での変数名は、キャメルケースでも無く、C語みたいにアンダーバーで結合するでもなく、ドットで結合するって方式です。

なんか今までの風潮だと、ドットはメソッドチェーンだとか、pythonみたいにパッケージ名と関数名の分離に使うと思ってた。慣れるとR語風が快適と思えるようになってきましたよ。

ここまで、たった4行書くだけで、上のようなデータが手に入り、R語の破壊力に感心しましたよ。これってベクトルを中心にするR語の威力なんですなぁ。

それから、emacsでESSを使ってる人用のTipです。実行したい行にカーソルを持って行って C-c C-cします。すると、空行に出会うまで、連続して実行してくれるんです。これで、途中々、実行結果を確認しながら、コードが書けます。何気ない気遣いに感謝んです。

ここまでが、前回からの復習。いよいよ時系列って事で、 Rと時系列(1) を、教科書にします。オブザーバーとして、実データを時系列解析して結果を考察

まずは、時系列形式を作るんだな。

hi.ts <- ts(m.am$hi, start=c(2012,1), frequency=12)
pl.ts <- ts(m.am$pl, start=c(2012,1), frequency=12)

データ列は、年月順に並べてあるもの。そしてデータの開始年、月を指定。月は12で繰り上がるよってヒントを与えてあげるんだな。

> str(hi.ts)
 Time-Series [1:104] from 2012 to 2021: 125 123 121 118 119 ...
> end(hi.ts)
[1] 2020    8

そのように解釈してくれたか確認。strって関数で要約を確認。普通strったら文字列関連かと思うけど構造ってのの略なんだな。で、年は2021年まで続いているって。構造の詳細を確認。そう、そこまでのデータが有るんよ。

> hi.ts
          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
2012 125.0645 123.2759 120.9667 118.1667 118.9032 115.8276 113.2903 117.5484 ..
2013 124.0323 119.7857 116.8065 117.0333 118.6154 120.6333 118.9355 119.4194
 :

一応、検算してみるかな。

 > colMeans( subset(m.am, yy == 2012 & mm == 8) )
         mm         yy         hi         lo         pl
    8.00000 2012.00000  117.54839   76.29032   58.90323
!> colMeans( subset(m.am, yy == 2013 & mm == 2) )
        mm        yy        hi        lo        pl
    2.0000 2013.0000  119.7857   75.0000   57.0000

どうやら問題なさそう。でも問題は、 こんな表を見せられたって、眼がチカチカするだけ。図にしましょ。

ts.plot(hi.ts, pl.ts)
plot(decompose(hi.ts))

最初の図は、取り合えずの出来具合を、同一の画内に表示。ふと思ったんだけど、最高血圧と脈拍って、単位が全く違う。そんな数値を同一に表示していいのか。

2番目の指示で、前回の教科書と同じ形式の図が出てきた。元データ、トレンド、周期性データ、誤差分。こうしてトレンドが年齢と共に上昇してるのが残酷なまでに表示されちゃったぞ。周期性では、1月が一番高く春が一番低いって結果。寒い時は気を付けましょうって事だな。

acf(hi.ts)
acf(hi.ts, lag.max=96,ylim=c(-.2,1.0))

自己相関。下のやつは、lagを広範囲に振った時。いずれも前月との強い相関が見られますた。

ccf(hi.ts, pl.ts)

自己相関だけじゃつまんないので、最高血圧と脈拍の相関も取ってみた。X軸はlagになってて、0のポイントで、強い負の相関って結果になった。まあ、上で並べて表示した時、予想はされていたけどね。

血圧が低いと脈拍を多くして、体内に送り出す血液量を増す必要がある。血管内にゴミが溜まって細くなると、圧力を高める必要が有る。血圧が上がると脈拍を低下させてエンジンである心臓をいたわろうと言う挙動になるのかな? 素人の推測だけど。

Nile

流量なんてのが出てきたんで、R語に備わっている流量データに注目してみる。有名なのはUKgasだ。これ、あちこちで使われているんで、少々秋田。で、他のを探したら、ナイル川の年間流量データってのが有った。

そんなに大きな川の流量なんて、どうやって測るねん? 年間って事だと、例えば1時間あたりのの流量を求めて、それを加算してくのかな。気が遠くなる。たとえ1時間の流量を求めるとしたら、どうやる? 根本原理を考察(素直に信じるのは良くないと思う)。

素人考えだけど、川の断面積を求める。川の流れるスピードを求める。それを掛け算すれば流量になるな。でも、流れは水面近くと底の方で違うだろうに。なんかそういうのの公式みたいなの有るのかな? 軽く検索したら、河川の土木事務所とかで使われているEXCELのアプリが有るみたい。公式なり、表を引いて計算してるんだろうね。

で、このナイル川の流量データには少々問題が有るようである。それを統計でカバーしましょって例みたい。

> example(Nile)

Nile> require(stats); require(graphics)

Nile> par(mfrow = c(2, 2))

Nile> plot(Nile)
Hit <Return> to see next plot:

Nile> acf(Nile)
 :

こんな風に、取り合えず実行してみてから、コードを眺めてみるのがいいのかな。これぞ現場実習です。で、現場ではついて行けず。なんせナイル川は遠すぎます。日本の河はどうよ。

ナイル川には負けるけど、日本で一番長い川、信濃川の流量データでもないかしら?

データで見る信濃川 って事で、信濃川スゲーってのが強調されてるけど、年間推移は見つからず。 ああ、ここに有った。日流量年表検索

なお世界一の河はアマゾン。日本もその流域の一部になってて、既に水没してますよ。中小企業は、軒並み被害にあってる。政府よ、アマゾン被害から企業を守る堤防を早く築けよ。

で、信濃川のデータは細かすぎる。年を指定してDLして、、、なんてやってられないね。まとめたやつは、日本河川協会からDVDになって発売されてる。昔から水の利権は争いの種だったけど、データにも利権が発生してるのね。

しょうがない河岸を変えよう。

日本の気温

同じ国土交通省の管轄だけど、気象庁はどうよ? 気象庁で分かり易いのは、気温だな。

よく天気予報のおねーさんが言ってる、先月は平年に比べて2度高かったとかの平年。 過去30年のデータを集めてきて、並べ替え。それを3等分して真ん中の組のデータだけを使って計算するとか。極端に寒い(暑い)データを除外した上で求めているのね。そして、このデータは、5年だか10年おきき更新するらしい。

こういう知識は気象予報士の試験に出るのかな。たまたま読んでいた小説に、リタイアした人が趣味で予報士試験に挑戦する話が出て来てた。難しくて3度の挑戦でやっと合格したとか。合格率5%とからしい。1000から1500時間の勉強が必要とか。暇を持て余した人はボケ防止に挑戦してみるのも面白い鴨。

ここにDL用のページが用意されてた。 過去の気象データ・ダウンロード

東京と富士山と言う、日本を代表する地域を指定。Fuji-yamaは、静岡県所属? それとも山梨県? よく出るクイズだな。観測所が有る所は山梨県でした。他に南極なんてのも有ったけど、無視。データの種類は、取り合えず温度データだけ。期間は1990年から先月までを指定。

結果を表示させて、問題なければCSVでダウンロード。優しいインターフェースで、小学生にも使えるぞ。夏休みの宿題のネタにどうぞって訳だな。

ちゃんとDL出来たか lireofficeで確認。文字コードをSHIFT-JISに指定しないと、文字化けするぞ。どこまで行ってもWindows仕様な訳ね。年・月を分離する記号は "/" になってるし、月も1桁だったり2桁だったり、EXCEL以外では不便な風。

ダウンロードした時刻:2020/09/10 05:30:27

,東京,東京,東京,富士山,富士山,富士山
年月,平均気温(℃),平均気温(℃),平均気温(℃),平均気温(℃),平均気温(℃),平均気温(℃)
,,品質情報,均質番号,,品質情報,均質番号
1990/1,5.0,8,1,-18.7,8,1
1990/2,7.8,8,1,-12.1,8,1
 :
2020/7,24.3,8,2,4.8,8,1
2020/8,29.1,8,2,7.6,8,1

生データを確認。品質情報とか付いてるけど、そんなのをオミットして取り込みたいぞ。 元データを改変する事なくね。R言語にはpipeが有る事を知った。

con <- pipe("sed '1,5d' data.csv | sed 's#/#,#' | cut -d',' -f1-3,6")
tmp <- read.csv(con, header=F)
colnames(tmp) <- c("yy", "mm", "tokyo", "fuji")

R言語でパイプから流れて来るデータを読めるのね。これは素敵な事。して、パイプでunix側に依頼したのは、まずsedで、余計な冒頭の5行を削除、データだけをむき出しにする。

次の段のsedで年・月を分ける"/"をカンマに変更。そうしておいて、cutコマンドで区切り記号をカンマに設定して、必要な項だけを抽出。

なお、2段目のsedで正規表現の区切りに"#"を使ったのは、見易くする為です。

!> summary(tmp)
        yy             mm             tokyo            fuji
  Min.   :1990   Min.   : 1.000   Min.   : 4.70   Min.   :-22.50
  1st Qu.:1997   1st Qu.: 3.000   1st Qu.: 9.20   1st Qu.:-14.55
  Median :2005   Median : 6.000   Median :16.90   Median : -5.70
  Mean   :2005   Mean   : 6.457   Mean   :16.53   Mean   : -5.94
  3rd Qu.:2012   3rd Qu.: 9.000   3rd Qu.:23.10   3rd Qu.:  2.45
  Max.   :2020   Max.   :12.000   Max.   :29.60   Max.   :  8.30
                                                  NA's   :1

東京暑いと言っても平均気温だと30度を割ってるのか。そんな事より気になるのは富士山のデータにNA(欠損値)が有る事。一体何時なんだ。

!> subset(tmp, is.na(tmp$fuji))
       yy mm tokyo fuji
 353 2019  5    20   NA

去年の5月か。全面的な改修工事でもあって、閉鎖されてたのかな?

sakae@pen:~/ANA$ grep '2019/5' data.csv
2019/5,20.0,8,2,,1,1

元データに当たると、確かにデータが欠損してた。品質情報(右から2番目)も、確かに、統計値がない(欠即)って説明されてたよ。

 > max(tmp$fuji)
 [1] NA
!> max(tmp$fuji, na.rm=T)
 [1] 8.3

欠損値が影響してくるような関数では、安全なために結果がNAになる。そのような関数には、NA除外オプションが用意されていたりする。NA値をどうするかってのは、悩ましい問題。色々な処理方法が提案されてる。例えば、Rによる欠損値処理 。やり出したらキリがないので先に進む。

tokyoのトレンドを表示させてみた。内心、じわじわと温度上昇してると思ったけど、結構長期 周期でふらふらしてた。

このトレンドって、観測データを平均化したもの(ですよね)。もう一回平均化したらどうなる?

 > x <- decompose(t.ts)$trend
!> plot(decompose(x))

t.tsはtokyoの時系列形式。これを一回処理してトレンドデータだけをxに取り出す。それをもう一回decompose。結果は、更に滑らかになった。って、当たり前の事だな。   今回はこれぐらいにするか。やや、新しいページを見つけたぞ。

Rではじめての時系列分析(AirPassengers)

このページを辿って行くと、pythonでは、どうやるなんてのに案内された。pythonめんどくせーな。

dplyr

検索してると、よく %>% と言う演算子に出会う。これは、とあるライブラリーとかで多用される パイプらしい。あるライブラリーの基本は、 dplyr — 高速data.frame処理 だったりするみたい。

ちょっとやってみる。

library(dplyr)
tmp %>%
    dplyr::group_by(mm) %>%
    dplyr::summarise(tokyo=mean(tokyo))

こんなやつを書いてみた。気温情報が信号源。そいつを月でグループ化する。その結果を集約させる。集約方法は、東京の平年気温を求めるって指定。中間結果が残らないので、綺麗。 結果は、

      mm tokyo
   <int> <dbl>
 1     1  6.08
 2     2  6.77
 3     3  9.86
 4     4 14.8
 5     5 19.3
 6     6 22.5
 7     7 26.3
 8     8 27.6
 9     9 24.1
10    10 18.8
11    11 13.5
12    12  8.63

スポットで、8月の値を求めて、答え合わせ

> colMeans( subset(tmp, mm == 8) )
         yy          mm       tokyo        fuji
2005.000000    8.000000   27.603226    6.351613

いやと言う程詳しい説明があった。 dplyr入門 (新版) これを使い込む程、のめり込む事ってあるのだろうか? 私たちのR: ベストプラクティスの探究 こういう、書きかけの本を見つけたぞ。価値ある逸品。

inside R

久しぶりにアプリの中身探求、って訳でdebian32Bit機で、Rをコンパイルしてみた。Verは3.5.2。debianにデフォで入っていたのと同じもの。

コンパイルに20分かかった。フォートランのコンパイラーとかが動いて大騒ぎですよ。入れた場所は、~/MINEの下。普通にPATHを通しておいて起動すると、~/MINI/bin/R と言うshellスクリプトが動いて、最終的には下記のようになる。Rのプロンプトが出てる時 help(filter)とかやった。自前のpagerが動いているな。

debian:~$ ps a | grep R
 2560 pts/2    S+     0:00 /home/sakae/MINE/lib/R/bin/exec/R
 3023 pts/2    S+     0:00 sh -c '/home/sakae/MINE/lib/R/bin/pager' < '/tmp/RtmpQKcM1v/a00ec65828'

動いている時、emacs-gdbを下記のように起動してアタッチした。

Run gdb (like this): gdb -i=mi /home/sakae/MINE/lib/R/bin/exec/R -p 2560

まずは、バックトレースして、どこら辺に居るか確認

(gdb) bt
#0  0xb7f5bd31 in __kernel_vsyscall ()
#1  0xb739f707 in waitpid () from /lib/i386-linux-gnu/libc.so.6
#2  0xb731d893 in ?? () from /lib/i386-linux-gnu/libc.so.6
#3  0xb74cfa08 in system_compat () from /lib/i386-linux-gnu/libpthread.so.0
#4  0x0062cd39 in R_system (command=0xbfd4ca10 "'/home/sakae/MINE/lib/R/bin/pager' < '/tmp/RtmpQKcM1v/a00ec65828'") at sysuti\
ls.c:330
#5  0x00683d4d in Rstd_ShowFiles (nfile=1, file=0x1c9c7e4, headers=0x1c9c808, wtitle=0x1ac5bc4 "R Help on \342\200\230filter\\
342\200\231", del=TRUE, pager=0x1c07e1c "/home/sakae/MINE/lib/R/bin/pager") at sys-std.c:1271
#6  0x005bd7bc in do_fileshow (call=0x22f7fe8, op=0xbc1d38, args=<optimized out>, rho=0x1cb6f68) at platform.c:365
#7  0x005641eb in bcEval (body=body@entry=0x2127eb8, rho=rho@entry=0x1cb6f68, useCache=useCache@entry=TRUE) at eval.c:6775
#8  0x005702db in Rf_eval (e=0x2127eb8, rho=0x1cb6f68) at eval.c:616
#9  0x00571e94 in R_execClosure (call=call@entry=0xb502fc4c, newrho=newrho@entry=0xbad9f8, sysparent=<optimized out>, rho=<op\
timized out>, arglist=<optimized out>, op=<optimized out>) at eval.c:1765
#10 0x00572efe in Rf_applyClosure (call=<optimized out>, op=<optimized out>, arglist=<optimized out>, rho=<optimized out>, su\
ppliedvars=<optimized out>) at eval.c:1693
#11 0x005685ab in bcEval (body=body@entry=0x14b2fd8, rho=rho@entry=0x2359738, useCache=useCache@entry=TRUE) at eval.c:6743
#12 0x005702db in Rf_eval (e=0x14b2fd8, rho=0x2359738) at eval.c:616
#13 0x00571e94 in R_execClosure (call=call@entry=0xb502f028, newrho=newrho@entry=0xbad9f8, sysparent=<optimized out>, rho=<op\
timized out>, arglist=<optimized out>, op=<optimized out>) at eval.c:1765
#14 0x00572efe in Rf_applyClosure (call=<optimized out>, op=<optimized out>, arglist=<optimized out>, rho=<optimized out>, su\
ppliedvars=<optimized out>) at eval.c:1693
(More stack frames follow...)
#15 0x005af446 in applyMethod (call=call@entry=0x2359798, op=op@entry=0x14b31f8, args=0x2359938, rho=0x23598d8, newvars=0x235\
97d8) at objects.c:118
#16 0x005b0030 in dispatchMethod (sxp=sxp@entry=0x14b31f8, dotClass=0x118f518, cptr=0xbfd4f578, method=0x17b0668, generic=0xc\
2313c "print", rho=0x23598d8, callrho=0x23599b8, defrho=0xbcc988, op=<optimized out>, op=<optimized out>) at objects.c:409
#17 0x005b04d2 in Rf_usemethod (generic=0xc2313c "print", obj=0x2368cac, call=0x14b26a8, args=0xbad9f8, rho=0x23598d8, callrh\
o=0x23599b8, defrho=0xbcc988, ans=0xbfd4e714) at objects.c:449
#18 0x005b07b4 in do_usemethod (call=<optimized out>, op=<optimized out>, args=<optimized out>, env=<optimized out>) at ../..\
/src/include/Rinlinedfuns.h:451
#19 0x0056435f in bcEval (body=body@entry=0x14b26c8, rho=rho@entry=0x23598d8, useCache=useCache@entry=TRUE) at eval.c:6159
#20 0x005702db in Rf_eval (e=0x14b26c8, rho=0x23598d8) at eval.c:616
#21 0x00571e94 in R_execClosure (call=call@entry=0xb502f010, newrho=newrho@entry=0xbad9f8, sysparent=<optimized out>, rho=<op\
timized out>, arglist=<optimized out>, op=<optimized out>) at eval.c:1765
#22 0x00572efe in Rf_applyClosure (call=<optimized out>, op=<optimized out>, arglist=<optimized out>, rho=<optimized out>, su\
ppliedvars=<optimized out>) at eval.c:1693
#23 0x005704a9 in Rf_eval (e=<optimized out>, rho=<optimized out>) at eval.c:739
#24 0x005d392e in Rf_PrintValueEnv (s=0x2368cac, env=<optimized out>) at print.c:1067
#25 0x005a1677 in Rf_ReplIteration (state=<optimized out>, browselevel=<optimized out>, savestack=<optimized out>, rho=<optim\
ized out>) at main.c:262
#26 Rf_ReplIteration (rho=<optimized out>, savestack=<optimized out>, browselevel=<optimized out>, state=<optimized out>) at \
main.c:198
#27 0x005a19bc in R_ReplConsole (rho=0xbcc9e8, savestack=0, browselevel=0) at main.c:308
#28 0x005a1a67 in run_Rmainloop () at main.c:1082
#29 0x004aa588 in main (ac=1, av=0xbfd509f4) at Rmain.c:29

Rのアプリが動くとネストが深くなる(深すぎる)。ってんで、下記はアイドル時

(gdb) bt
#0  0xb7f5bd31 in __kernel_vsyscall ()
#1  0xb73d1d21 in select () from /lib/i386-linux-gnu/libc.so.6
#2  0x006827e9 in R_SelectEx (n=1, readfds=0x8f2840 <readMask.14050>, writefds=0x0, exceptfds=0x0, timeout=0x0, intr=0x6825e0\
 <handleInterrupt>) at sys-std.c:150
#3  0x00682af3 in R_checkActivityEx (usec=-1, ignore_stdin=0, intr=0x6825e0 <handleInterrupt>) at sys-std.c:327
#4  0x00683092 in Rstd_ReadConsole (prompt=0xbc9f10 "> ", buf=0xbfd4f8e8 "\n", len=4096, addtohistory=1) at sys-std.c:1005
#5  0x005a1427 in Rf_ReplIteration (rho=0xbcc9e8, savestack=0, browselevel=0, state=0xbfd4f8dc) at main.c:206
#6  0x005a19bc in R_ReplConsole (rho=0xbcc9e8, savestack=0, browselevel=0) at main.c:308
#7  0x005a1a67 in run_Rmainloop () at main.c:1082
#8  0x004aa588 in main (ac=1, av=0xbfd509f4) at Rmain.c:29

この2つを見比べながら、潜って行けばいいんだな。

R Internals を見てくれか、な。英語は読書スピードが落ちるんで翻訳しながら読んで行こう。


This year's Index

Home