R言語 (3)

from source
- startup
- main
- datasets
- tcltk
- po
欠損値の処理
暑さ寒さは彼岸まで？
input

from source

前回はsourceからRを作って、恒例のgdbにかけてみた。それで終わってはもったいないので、少し寄り道する。

startup

libの中にご本尊のRが鎮座してるので、それを直接起動

debian:~$ MINE/lib/R/bin/exec/R
MINE/lib/R/bin/exec/R: error while loading shared libraries: libRblas.so: cannot open shared object file: No such file or directory

ライブラリーの一部が見つからないと文句を言われた。ならば伝統のライブラリィーパスの指定だな。

debian:~$ ls MINE/lib/R/lib
libRblas.so*  libRlapack.so*
debian:~$ export LD_LIBRARY_PATH=/home/sakae/MINE/lib/R/lib
debian:~$ MINE/lib/R/bin/exec/R
Fatal error: R home directory is not defined

少しは先に進んだけど、今度は別のエラーが発生

debian:~$ export R_HOME_DIR=/home/sakae/MINE/lib/R
debian:~$ export R_HOME=$R_HOME_DIR
debian:~$ MINE/lib/R/bin/exec/R

R version 3.5.2 (2018-12-20) -- "Eggshell Igloo"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: i686-pc-linux-gnu (32-bit)
 :
Error in library(ggplot2) : there is no package called ‘ggplot2’
[Previously saved workspace restored]
>

bin/Rって起動スクリプトを参考に、必要そうな環境変数を設定。これでやっと起動してきた。 ggplot2が見つからないエラーは、.Rprofileで指定されたのが参照されたんだな。でもアドオンの場所に、そんなの無かったと言ってきてるんだ。こういう遊びをやっておくと、後で問題が発生してもオロオロする事が無いと思うぞ。

main

Rを構成している本体は、src/mainの中に鎮座してる。ほとんどがC語なんだけど、一つだけxxxpr.fってフォートラン語が混じっていた。ここの場所で、

debian:main$ etags *.[ch]

しとけば、軽やかにソース内を飛び割れるぞ。

datasets

Rの本体以外は、libraryってdirの中に格納されてる。datasetsって言う楽しい名前のやつも有った。データもRの一部ですってのが嬉しい。興味深く覗いてみるか。

dataとmanの部に分かれていた。その中のdataでは、何時もお世話になる　iris.Rが居た。統計の父フィッシャー先生の研究データ。

"iris" <- data.frame(
Sepal.Length = c(5.1, 4.9, 4.7, 4.6, 5, 5.4, 4.6, 5, 4.4, ...),
Sepal.Width = c(3.5, 3, 3.2, 3.1, 3.6, 3.9, 3.4, 3.4, 2.9, ...),
Petal.Length = c(1.4, 1.4, 1.3, 1.5, 1.4, 1.7, 1.4, 1.5, 1.4, 1.5, 1.5, ...),
Petal.Width = c(0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.3, 0.2, ...),
Species = gl(3,50, label = c("setosa", "versicolor", "virginica")))

こんな風にデータを持ってる訳ね。次は、代表的な時系列データを見ておく。

Nile.R

Nile <- stats::ts(c(1120, 1160, 963, 1210, 1160, 1160, 813, ...),
start=1871)

UKgas.R

UKgas <- stats::ts(c(160.1, 129.7, 84.8, 120.1, 160.1, ...),
start = 1960, frequency = 4)

最後はおまけで、R発案者の母国ニュージーランドの国土地理院発行？の、volcano火山のデジタル地図（標高）を見る。

"volcano" <-
structure(c(100, 101, 102, 103, 104, 105, 105, 106, 107, ...),
.Dim = c(87, 61))

Rって地理系もサポートしてるんですね。面白そう。

tcltk

同じくライブラリーの一員としてtcltkが有った。グラフィックが今のように発達していなかった頃には、重宝されました。昔のよしみで正会員の特権を保持してます。

demo(tkfaq)

とかして、眺めるもよし、下記の簡易GUIをやってみるのも、昭和の時代にタイムスリップしたようで楽しいかも。

library(tcltk)
tt <- tktoplevel()
tkpack(txt.w <- tktext(tt))
tkinsert(txt.w, "0.0", "plot(1:10)")

eval.txt <- function()
    eval(parse(text = tclvalue(tkget(txt.w, "0.0", "end"))))
tkpack(but.w <- tkbutton(tt, text = "Submit", command = eval.txt))

テキストエリアにplotコマンド（のテキスト）を表示しておいて、送信ボタンで、それを実行する例。使い終わった、 tkdestroy(tt)　して、痕跡を消しておく事。

pb <- tkProgressBar("test progress bar", "Some information in %",
                    0, 100, 50)
Sys.sleep(0.5)
u <- c(0, sort(runif(20, 0, 100)), 100)
for(i in u) {
    Sys.sleep(0.1)
    info <- sprintf("%d%% done", round(i))
    setTkProgressBar(pb, i, sprintf("test (%s)", info), info)
}
Sys.sleep(5)
close(pb)

プログレスバーの表示方法。CUIの人も、これぐらいなら使ってもいいかも。

po

ソースファイルの要所にpoって名前の翻訳済み文字列ファイルが用意されてる。我が日本もある。どなたかの貢献であろう。普段は英語版で使っているんだけど、たまには、日本語してみるかな。

debian:tmp$ LANG=ja_JP.utf8 R
 :
'demo()' と入力すればデモをみることができます。
'help()' とすればオンラインヘルプが出ます。
'help.start()' で HTML ブラウザによるヘルプがみられます。
'q()' と入力すれば R を終了します。

欠損値の処理

前回の富士山温度データに欠損が有った。

352 2019  4  13.6 -14.9
353 2019  5  20.0    NA
354 2019  6  21.8   5.3
355 2019  7  24.1   5.2
356 2019  8  28.4   6.9
357 2019  9  25.1   5.1

欠損値は無視するって言う na.rm=T が有ったけど、こんなの手軽に使っていいのか？　今回の場合は、問題の５月のデータは正規で30個。それが一つ欠けても大勢に影響は無いだろう。統計的な裏付けとして、30個の母数から29個をサンプリングして平均を求めましたって考え。

気味が悪いので、欠けたデータをねつ造する方法を考える。すぐに思い付くのは、4,6月の平均を使っちゃう方法。一番簡単そう。

でも、欠けたデータが８月だとしたらどうよ？　８月は一年で一番気温が上昇する月。7,9月の平均を取っちゃうと、2度も冷夏な年になっちゃう。さて、どうする？

思い付くのは、29年間のデータの平均を求めて、それを埋め込んじゃ方法。

!> mean(subset(tmp$fuji, tmp$mm == 8))
 [1] 6.351613

欠損値は余り気にする必要は無いってスタンスが多いけど、場合によっては大変になる事があるそうな。そんな時の解決方法が、欠損データの処理で、解説されてた。

暑さ寒さは彼岸まで？

この所、涼しくなったかと思うと、また暑さがぶり返したりしていて、秋と夏が喧嘩してるみたい。昔の人は、そんな状況を、「暑さ寒さは彼岸まで」と言っていた。本当か？

「暑さ寒さも彼岸まで」は本当？　データで検証してみた

二番煎じの感はあるけど、オイラーもやってみる。

気象庁へ行って、ローカルエリアのデータを取ってきた。（URLは前回参照）

2010-2019年にかけての9月と10月の日毎の平均、最高、最低温度。これだけ有れば何とかなるだろう。

library(dplyr)

con <- pipe("sed '1,6d' data.csv | sed 's#/#,#g' | cut -d',' -f1-4,7,10")
so  <- read.csv(con, header=F)
colnames(so) <- c("yy", "mm", "dd", "av", "hi", "lo")

取ってくるデータによってCSVの項目位置が微妙に違うので、前回からちょびっと修正した。一つ注意が有る。EXCEL仕様の 2010/9/20 みたいなのの "/" を "," に置き換える時、正規表現の最後に g を付けておく事だ。このgが無いと、最初の / を変換しただけで満足してしまう。

取ってきたデータに欠損値が無いか一応確認。普通はsummary(so)とか使うんだけど、それだと、yy,mm,ddも出てきちゃう。そんな訳で、必要なものだけを選んでサマリーしてみた。 Rが提供するpipe関数といい、libraryがサポートしてるパイプといい、なかなか便利。

> so %>% select(av,hi,lo) %>% summary()
        av              hi              lo
  Min.   : 5.50   Min.   : 7.70   Min.   : 1.80
  1st Qu.:14.30   1st Qu.:19.50   1st Qu.:10.20
  Median :18.20   Median :23.20   Median :14.15
  Mean   :17.99   Mean   :23.28   Mean   :13.95
  3rd Qu.:21.40   3rd Qu.:26.98   3rd Qu.:17.80
  Max.   :28.50   Max.   :34.80   Max.   :24.90

わずか２月の間に、気温が大幅に低下してる。まあ、こちらには負けるけどね。米コロラド州で48時間の間に気温38.3℃から1.1℃に急変。

冒頭で利用するライブラリーを宣言するのは、良い習慣だ。で、今回のやつを宣言すると

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

こんな警告が出て来る。名前の衝突が発見されたんで、stats::{filter,log}とかは、マスクしとな。statsとかbaseにあるやつを使いたい場合、パッケージ名を前置すればいいんだな。

こういう事が有るから、最初にロードしとくって重要。

データは時系列（9,10月が年毎に繰り返されるるやつ）なんで、それ風にしとく。

av.ts <- ts(so$av, start=c(2010,9),frequency=61)
hi.ts <- ts(so$hi, start=c(2010,9),frequency=61)
lo.ts <- ts(so$lo, start=c(2010,9),frequency=61)

61日分のデータが続いたら、年がインクリメントされるって指定で大丈夫だったよ。

ts.plot(window(av.ts,c(2012,1),c(2012,61)))

普通にグラフを書いちゃうと10年分が表示されちゃう。そういう時はwindowを使って期間を抜き出せばよい。但し、X軸の表示がチト怪しいけどね。

上で時系列データを作ったけど、こういう場合は使い道無いな。素直にウェザーニューズに倣ってやってみるか。それも能の無い話だ。

数字で見せられるよりグラフでしょう。統計っぽく、Rによる箱ひげ図の描画を参考にね。

show <- function() {
    for (i in unique(so$yy) ){
        q <- so %>% filter(yy== i) %>% select(av)
        boxplot(q[10:19,], q[20:29,], q[30:39,], main = as.character(i))
        dmy <- readline("Hit <RET> key to next ")
    }
}

元記事では、彼岸の前後の週も含めて、平均値が列挙されて表で結果が示されている。でも平均値と言うと丸められてしまって、寒暖差が見えない。ならば、箱ひげ図だろう。

一つの箱ひげデータが１週間では、ちと心もとないと思ったので、10日区切りにした。一年毎に１グラフ。年を抽出する為、ユニーク関数を使ってみた。年でフィルターして、平均気温の61個のデータをｑに保存。

行番号は連番になってるんで、それをインデックスにして、抜き出した。最初、中間変数ｑを使うのを止めて、パイプで繋ぎ、boxplot( .[10:19,], .[20:29,], ..) としたけど、最初の箱しかグラフにならなかった。パイプで来たデータは展開されないのね。なんかそこを逃げる方法が有りそうだけど、深くは追及しない。

boxplotのタイトルは、main って引数に文字列で指定するようだ。helpでの説明には出ていなかったぞ。例を実行したら、使用してる奴が居て、それを真似させてもらった。それから、数字を文字列にする方法で、ちょっと悩んだのは秘密だ。

mainは何者ってんで、boxplot.Rを参照してみると、こんなコードが有った。

do.call("title",
        pars[names(pars) %in% c("main", "cex.main", "col.main",
                                "sub", "cex.sub", "col.sub",
                                "xlab", "ylab", "cex.lab", "col.lab")])

グラフィック関係の山ほど有るオプションの一員ぽい。全くGUI系は難儀だなあ。

まあ、これで、年毎にグラフを生成させて、フムフム出来るようになったけど、2019年のデータを見てて、あれ？　前年はどうだったって希望は叶えてくれない（老人性痴ほう症）。

ならば一遍に見られるようにしちゃえ。

# arg k for what graph: "lo", "av" (default), "hi"
mkgr <- function(k="av") {
    par(mfrow = c(2,5))
    for (Y in unique(so$yy) ){
        q <- so %>% filter(yy== Y) %>% select(all_of(k))
        boxplot(q[10:19,], q[20:29,], q[30:39,],
                main = as.character(Y), ylim = c(5,35))
    }
}

mfeowってパラメータを宣言して、縦に２段、横に５つのグラフ領域を指定。この描画デバイスにグラフを書くと、勝手にエリアを埋めてくれる。またboxplot中にylimを設定して、表示温度範囲を固定した。こうしておかないと、それぞれのグラフ内で最適化されてしまって、比較にならない。

この温度範囲の決定にあたって、最初は c(10,30) にしてたんだけど、これだと、最高気温、もしくは最低気温でレンジ外になる年が有ったんで、広めた。こんな所でも、この時期は寒暖の差が激しいって分かって、面白かったぞ。

> mkgr(k='hi')
 Note: Using an external vector in selections is ambiguous.
  Use `all_of(k)` instead of `k` to silence this message.
  See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
 This message is displayed once per session.

それから、最初debian機のR3.5.2で、パイプラインの最後のselectのを、select(k)ってしてたんだ。たまにはウブでもと思ってR3.6.4でも試したら、上記のような教育的指導がなされた。URL先にはFAQって書いてあるんで、みんなこの道を辿っているのね。安心した。

今度は別の要求。最高気温、平均気温、最低気温を同時に眺めてみたい。これはもうPDFにでも落として、印字してから眺める鹿。

pdf()
mkgr(k='hi')
dev.off()

こうすると、Rplots.pdfって名前のファイルが出来上がる。何時も同じ名前じゃいやだって人はpdfデバイスを開く時に、ファイル名を指定してもよい。更に、紙の大きさとかfontsの指定等も、pdf関数の中で指定出来る。

これをUSBに焼いて、散歩の途中でコンビニに寄ってくれば良い。

Rplots.pdf

改めて説明するまでも無いと思うけど、それぞれの年度内では、左から、９月中旬、９月下旬（アバウトに彼岸期間中）、１０月上旬って並びになってる。

箱ひげの見方。データを昇順に並べる。上下の横棒は、最高・最低値。箱の上辺はデータの個数の3/4ポイントの値、下辺は1/4ポイントの値。箱の中の太い横棒は、丁度真ん中の値（ほぼ平均値）となる。横棒の外側に有る丸は、外れ値（平均から大幅にずれたやつ）を表している。

拡がりが大きいほど、寒暖差がある、体に堪えるって事だ。

今のコードは、落としてきたデータの半分しか使っていない。全部のデータを使ってやらないともったいない。各月を前半・後半に分けて、４本の箱ひげセットにしてみるかな。ひと箱、15日分のデータになるんで、箱の威力を増すだろう。改造は直ぐに出来るからやってみるか。そして、春の彼岸も確認してみるかな。

とか思っていたら、女房が言うには、この時期（春の彼岸も）、頭が痛くなるそうな。気温より気圧が原因と言うけど、どうだろう。気圧も含めて検証しろって事かな？

夏の異常気象をオープン・データで確認 (1/2) あっ、こんなのも有った。

input

上でグラフを一枚ずつ送って行く方法。多分pauseとかだろうってんで、血眼になって探したけど、見つからず。そんな検索の途中で見つけたやつ。

> q <- menu(c('R', 'scheme', 'python'))

1: R
2: scheme
3: python

Selection: 2
> q
[1] 2

なんかCUIのI/Fで良く使ってそうなのが有った。脳の肥やしなるな。

で、件のやつは、素直にreadline()で良いらしい。

でも、グラフを書く exampleを実行すると、 "Hit <Return> to see next plot:"　こんなのが出てくる。何処で出しているのだろう？　ちょっと気になるんで探してみる。

./library/base/po/ja.po:msgid "Hit <Return> to see next plot: "
./library/base/po/nn.po:msgid "Hit <Return> to see next plot: "
  :
./main/devices.c:       R_ReadConsole(_("Hit <Return> to see next plot: "), buf,
 1024, 0);

poファイルは現地語用の翻訳データが入ったもの。原本は英語で書かれている。 main/devices.c

void NewFrameConfirm(pDevDesc dd)
{
    if(!R_Interactive) return;
    /* dd->newFrameConfirm(dd) will either handle this, or return
       FALSE to ask the engine to do so. */
    if(dd->newFrameConfirm && dd->newFrameConfirm(dd)) ;
    else {
        unsigned char buf[1024];
        R_ReadConsole(_("Hit <Return> to see next plot: "), buf, 1024, 0);
    }
}

この関数、何処から呼ばれてる？　追ってみたら、library/graphics/src/graphics.c からだった。更に潜ったら面白いかも。

This year's Index

Home