R言語 (4)
外れ値
気象庁から頂いてきた、彼岸あたりの温度変化（3,4月 or 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")
mkgr <- function(k="av") {
par(mfrow = c(2,5))
for (Y in unique(so$yy) ){
q <- so %>% filter(yy== Y) %>% select(k)
boxplot(q[1:15,], q[16:30,], q[31:45,], q[46:60,],
main = as.character(Y), ylim = c(-10,30))
}
}
for (K in c('hi', 'av', 'lo')) {
pdf( paste(K, '.pdf', sep='') )
mkgr(k=K)
dev.off()
}
出来上がったPDFファイルをブラウザーで全て読み込んでおく。そうしておいて、Ctl+TAB で、画面を素早く切り替えれば、パラパラ漫画の出来上がり。わざわざコンビニで印刷して来る事は不要で、楽チンである。
春の彼岸を見ているんだけど、箱ひげ図の上下バーを外れて外側にある値が、小さい丸で表現される、外れ値がやはり有るなあ。並みの言葉で言えば、気温の特異値と解釈出来る。
箱ひげ図（Rが書いてくれるんだけど）の効用として、"データの統計分布での外れ値を発見すること、分布が不均質ではないかを確認することなどのための視覚化表現"って某Webに説明されてた。まさにその通りだな。
良いサンプルが出ていた。
set.seed(101) z <- rnorm(100000) boxplot(z) hinges <- qnorm(c(0.25,0.75)) IQR <- diff(qnorm(c(0.25,0.75))) abline(h=hinges,lty=2,col=4) ## hinges ~ quartiles, blue line abline(h=hinges+c(-1,1)*1.5*IQR,col=2) # red line
やってる事が良く分かった（積り）。で、詳細はソース嫁。
> boxplot
function (x, ...)
UseMethod("boxplot")
<bytecode: 0x55d34716e808>
<environment: namespace:graphics>
んが、出てこない。こういう時は、boxplot. に続いてTABを叩き、関係者を教えて貰え。
> boxplot.stats
function (x, coef = 1.5, do.conf = TRUE, do.out = TRUE)
{
if (coef < 0)
stop("'coef' must not be negative")
nna <- !is.na(x)
n <- sum(nna)
stats <- stats::fivenum(x, na.rm = TRUE)
iqr <- diff(stats[c(2, 4)])
if (coef == 0)
do.out <- FALSE
else {
out <- if (!is.na(iqr)) {
x < (stats[2L] - coef * iqr) | x > (stats[4L] + coef *
iqr)
}
else !is.finite(x)
if (any(out[nna], na.rm = TRUE))
stats[c(1, 5)] <- range(x[!out], na.rm = TRUE)
}
conf <- if (do.conf)
stats[3L] + c(-1.58, 1.58) * iqr/sqrt(n)
list(stats = stats, n = n, conf = conf, out = if (do.out) x[out &
nna] else numeric())
}
<bytecode: 0x563d0b91ace0>
<environment: namespace:grDevices>
新 「 数 学 I 」 統 計 分 野 に 於 け る 「 Ｒ 言 語 」 支 援 に よ る 教 材 作 成 の 一 例 と し て
学校でもRを使って箱ひげ図を勉強するのか。パソコンと言ったらEXCELって風潮に風穴を開けてもらいたいぞと。
EXCELはユーザーが付きっきりになって操作をせにゃならん。上で挙げたRのスクリプトと同じ事をしようとすれば、切った貼ったを幾千回してやっと、箱ひげ用データが揃う。違うデータを分析しようと思ったら、また同じ事の繰り返し。EXCELは米国防省だか商務省が送り込んだ、日本の生産性を低下させる兵器さ。
まて、EXCELでの操作を録画して、後で繰り返す事だって出来るぞ。他にもマクロを書くって方法だって有るぞ。まあ、好きにしてくださいな。
debugger
ってな訳でRを推すんよ。曲がりなりにもRは言語だ。虫が湧いた時の事を事前に調べておｋ。ユーザーレベルのデバッグ方法について調べてみた
ふむ、エラーになったら取り合えず、 traceback() させればいいんだな。そして虫の潜む関数が判明したら debug(hoge) しろとな。
早速、冒頭のやつでやってみる。最後のforの所を、main()って関数に収容しておいたからね。それとBugを仕込んでおいた。
> main() Error in xj[i] : only 0's may be mixed with negative subscripts
エラーで止まったのは分かるけど、ちと不親切。
!> traceback()
6: `[.data.frame`(q, -16:30, ) at ana.R!uHAHRE#5
5: q[-16:30, ] at ana.R!uHAHRE#5
4: boxplot.default(q[1:15, ], q[-16:30, ], q[31:45, ], q[46:60,
], main = as.character(Y), ylim = c(-10, 30))
3: boxplot(q[1:15, ], q[-16:30, ], q[31:45, ], q[46:60, ], main = as.character (Y),
ylim = c(-10, 30)) at ana.R!uHAHRE#5
2: mkgr(k = K) at ana.R!zb8S0f#4
1: main()
これで、何かあると何時も不機嫌に文句を垂れてくるpython並みになったな。もう虫の所在は分かってしまったけれど、一応ストーリー仕立てなんで。
> debug(boxplot.default)
!> main()
debugging in: boxplot.default(q[1:15, ], q[-16:30, ], q[31:45, ], q[46:60,
], main = as.character(Y), ylim = c(-10, 30))
debug: {
args <- list(x, ...)
namedargs <- if (!is.null(attributes(args)$names))
:
}
Browse[2]> c
Error in xj[i] : only 0's may be mixed with negative subscripts
>
引数で渡ってきたvecterの範囲がマイナスになってますねぇ（わざとらしいエラー）。まあ、作り付けの関数内にBugなんて事は無いでしょうから、自分で書いたコードを眼を皿にして追う。散歩に出かけてリフレッシュしてくるが、Bug取りの最高の手段だと思うよ。
debug関数は、ターゲットの関数の冒頭にbrowse()を埋め込むんだな。だからbrowse関数のコマンドが使える。n,c,Q ぐらいしか使い道が無いけどね。偽Bugを潰して、boxplotの挙動をステップ実行していける。変数も確認出来るから便利
Browse[2]> n debug: if (is.null(attr(groups, "names"))) attr(groups, "names") <- 1L:n Browse[2]> debug: names <- attr(groups, "names") Browse[2]> n debug: cls <- vapply(groups, function(x) class(x)[1L], "") Browse[2]> names [1] "" "" "" "" !Browse[2]> x [1] 34.1 34.3 32.9 32.2 34.8 34.2 32.9 24.4 30.0 30.3 31.9 28.3 29.5 28.3 25. 9
使い終わったら、undebug(hoge) しておこう。debugonce(hoge)って言う、一度だけdebuggerを呼び出すのも便利だぞ。
範囲
今回初めて気が付いた。いわゆる範囲問題。
> x <- 1:5 > x [1] 1 2 3 4 5 > x[0:7] [1] 1 2 3 4 5 NA NA > x[-1:4] Error in x[-1:4] : only 0's may be mixed with negative subscripts
タイトルを、範囲なんて事で始めちゃって、範囲にはマイナスを取れないと、一人合点してんだけど、どうやら、超赤恥ものだった。
マイナスは、指定場所の削除を意味するのね。【R言語】行列事始
> x <- 1:5 > x[-3] [1] 1 2 4 5
答え一発
ちょいと打ちひしがれて、別の道に走る。上の方でEXCELを貶してしまったけど、Rはどうよ？ Rを起動して、source(…) して、q() じゃ、まどろっこしい。
答え一発、コマンドラインからの入力だけで仕事を終えたい。これが出来れば、EXCEL連中の鼻をへし折れるってもんです。
sakae@ub:~/ANA$ Rscript ana.R
Attaching package: ‘dplyr’
:
yy mm dd av hi lo
609 2019 10 30 13.3 19.0 10.5
610 2019 10 31 12.8 19.4 8.3
スクリプトの中へ tail() を忍ばせて、ちゃんと動いているか確認させた。 パッケージのロードで、何たらかんたら言ってくるのがうざいな。
sakae@ub:~/ANA$ Rscript ana.R >/dev/null 2>&1
これだと、全く静かにさせる事が出来るな。でも、確認用も闇に葬られた。ライブラリーのメッセージだけ消したい。helpしろ。答えは
library(dplyr, warn.conflicts=F)
これで無様なメッセージは抑止出来る。でも臭い物には蓋って態度でいいのか？
元々有った機能をマスクしてまで、ユーザーが指定した関数（ライブラリーのロードだけど）が選ばれているんだから、作者がさぞ気合を入れて書いたものなんでしょう。
debugの有効活用
どんな風に気合が入ってるかソースを眺めてみたい。で、 debug(filter)とかの登場になる。
> source('ana.R')
> debug(filter)
> mkgr()
debugging in: filter(., yy == Y)
debug: {
UseMethod("filter")
}
Browse[2]> n
debug: UseMethod("filter")
Browse[2]>
debugging in: filter.data.frame(., yy == Y)
debug: {
as.data.frame(filter(tbl_df(.data), ..., .preserve = .preserve))
}
何か、日が暮れそうな雰囲気だな。
CRAN – search 経由で、 dplyr: A Grammar of Data Manipulation にたどり着く。後はtar玉を落としてきて、、、filter.Rを拝見
filter <- function(.data, ..., .preserve = FALSE) {
UseMethod("filter")
}
#' @export
filter.data.frame <- function(.data, ..., .preserve = FALSE) {
loc <- filter_rows(.data, ...)
dplyr_row_slice(.data, loc, preserve = .preserve)
}
こんな、書式のコメントが有った。helpだかhtmlの原稿になるのかな？
#' @family single table verbs
#' @inheritParams arrange
#' @param ... <[`data-masking`][dplyr_data_masking]> Expressions that return a
#' logical value, and are defined in terms of the variables in `.data`.
#' If multiple expressions are included, they are combined with the `&` operator.
#' Only rows for which all conditions evaluate to `TRUE` are kept.
#' @param .preserve Relevant when the `.data` input is grouped.
#' If `.preserve = FALSE` (the default), the grouping structure
#' is recalculated based on the resulting data, otherwise the grouping is kept as is.
#' @return
#' An object of the same type as `.data`. The output has the following properties:
:
pythonみたいに、関数の中に、 ''' hoge ''' みたいになってなくて、コメント酔い防止と言う意味では（オイラー敵に）好ましいと思うぞ。
つらつら見てたら
> starwars # A tibble: 87 x 13 name height mass hair_color skin_color eye_color birth_year gender <chr> <int> <dbl> <chr> <chr> <chr> <dbl> <chr> 1 Luke… 172 77 blond fair blue 19 male 2 C-3PO 167 75 NA gold yellow 112 NA 3 R2-D2 96 32 NA white, bl… red 33 NA 4 Dart… 202 136 none white yellow 41.9 male 5 Leia… 150 49 brown light brown 19 female 6 Owen… 178 120 brown, gr… light blue 52 male 7 Beru… 165 75 brown light blue 47 female 8 R5-D4 97 32 NA white, red red NA NA 9 Bigg… 183 84 black light brown 24 male 10 Obi-… 182 77 auburn, w… fair blue-gray 57 male # … with 77 more rows, and 5 more variables: homeworld <chr>, species <chr>, # films <list>, vehicles <list>, starships <list>
こんな楽しいデータを提供してた。
それから、このパッケージはWindowsユーザーに手厚いサポートをしてるように見受けられる。 まあ、EXCEL代わりにR言語を使いましょってのが主たる理由でしょ。そうすると、Rstadioを標準で選ぶだろう。そういうGUIへの親和性は大切ってのは、頷ける所です。
paste0
冒頭の例で、文字列の結合に paste を使った。このpasteには簡易版のpaste0が用意されている。何が簡易かと言うと、helpしろでもいいんだけどソース嫁で行ってみる。
debian(64Bit版）に、R4.0.2を取ってきて入れたんで、それをみる。configure時にX11が無いだとかpcore2が無いだとか抜かすんで、そんなの無いモードにした。make -j4 の威力は強大で、5分少々でコンパイル完了。32BitのシングルCPUでは30分かかったからね。
library/base/paste.R
## change in R 4.0.1: do it this way to allow for calls without new arg
## to work from earlier versions of R.
paste <- function (..., sep = " ", collapse = NULL, recycle0 = FALSE)
{
if(isTRUE(recycle0))
.Internal(paste(list(...), sep, collapse, recycle0))
else
.Internal(paste(list(...), sep, collapse))
}
paste0 <- function(..., collapse = NULL, recycle0 = FALSE)
{
if(isTRUE(recycle0))
.Internal(paste0(list(...), collapse, recycle0))
else
.Internal(paste0(list(...), collapse))
}
細かい修正が入っているんですなあ。問題はそこではなくて、C語を呼び出している点。
こんなのを実行
> paste0('hoge', 'fuga')
[1] "hogefuga"
gdb側ではBPを張っておくも、アドレス未解決。当然ヒットせず。
(gdb) info breakpoints Num Type Disp Enb Address What 1 breakpoint keep y <PENDING> paste0 2 breakpoint keep y <PENDING> _Internal
btすると
#3 0x000055f14bd46584 in Rstd_ReadConsole (prompt=<optimized out>, buf=0x7ffc2\
6de1aac "paste0('hoge', 'fuga')\n", len=4096, addtohistory=<optimized out>) at \
sys-std.c:1039
こんな具合に、入力したのが残っている。はてな？である。
んな訳で、 Rstd_ReadConsole にBPを置いてから違う文字列を評価すると、
> paste0('foo', 'bar')
[1] "foobar"
Rの画面には、先に答えが出ている。
そしてgdbの方は、
Breakpoint 3, Rstd_ReadConsole (prompt=0x55f14bfcb7f0 "> ", buf=0x7ffc26de1aac \
"paste0('foo', 'bar')\n", len=4096, addtohistory=1) at sys-std.c:952
952 if(!R_Interactive) {
どうも、libreadlineが途中で介入してて、かつ割り込みハンドラーも関係してるようで、素直に動かんなあ(一見不思議な現象は、多分readlineが残骸を残しているのだろう）。
Rf_ReplIteration にBPを置いて、地道に追って行くかな。
翌朝、たわわに実った田んぼを散歩してたら、神の声が聞こえた。paste0にBPを貼った時pendingって言われたろう？ 何故か考えるのじゃと。
考えたら、mein側にそんな関数はロードされていない、そもそもそんな関数は無いからgdbは対処しようがありませんって事だな。
ならば、C語のエリアでpaste0とかを探してみよう。そして行き着いた所は、paste.c。冒頭にこうある。
/* .Internal(paste (args, sep, collapse)) .Internal(paste0(args, collapse)) * do_paste uses two passes to paste the arguments (in CAR(args)) together. * The first pass calculates the width of the paste buffer, * then it is alloc-ed and the second pass stuffs the information in. */
最初に引数である文字列の長さを確認。次に長さに見合う場所を確保して、連結するとな。
早速、 do_paste にBPを置いて実行。
(gdb) bt #0 do_paste (call=0x9f2c48, op=0x928888, args=0x17cf3d8, env=0x17cf458) at paste.c:53 #1 0x005441eb in bcEval (body=body@entry=0x9f0c88, rho=rho@entry=0x17cf458, useCache=useCache@entry=TRUE) at eval.c:6775 #2 0x005502db in Rf_eval (e=0x9f0c88, rho=0x17cf458) at eval.c:616 #3 0x00551e94 in R_execClosure (call=call@entry=0xb4fb6010, newrho=newrho@entry=0x91b9f8, sysparent=<optimized out>, rho=<optimized out>, arglist=<optimized out>, op=<optimized out>) at eval.c:1765 #4 0x00552efe in Rf_applyClosure (call=<optimized out>, op=<optimized out>, arglist=<optimized out>, rho=<optimized out>, suppliedvars=<optimized out>) at eval.c:1693 #5 0x005504a9 in Rf_eval (e=<optimized out>, rho=<optimized out>) at eval.c:739 #6 0x005815ee in Rf_ReplIteration (state=<optimized out>, browselevel=<optimized out>, savestack=<optimized out>, rho=<optimized out>) at main.c:258 #7 Rf_ReplIteration (rho=<optimized out>, savestack=<optimized out>, browselevel=<optimized out>, state=<optimized out>) at main.c:198 #8 0x005819bc in R_ReplConsole (rho=0x93a9e8, savestack=0, browselevel=0) at main.c:308 #9 0x00581a67 in run_Rmainloop () at main.c:1082 #10 0x0048a588 in main (ac=1, av=0xbfeda254) at Rmain.c:29
おお、捕まえたぞ（ライ麦畑で捕まえてじゃなくて、田んぼの神様の思し召し。なんたってここは、お米の国ですから）
if(use_sep) { /* paste(..., sep, .) */
sep = CADR(args);
:
} else { /* paste0(..., .) */
=> u_sepw = sepw = 0; sep = R_NilValue;/* -Wall */
collapse = CADR(args);
}
これだけ見せられたら、lispのインタープリターだな。安心したぞ。
フレームを少々遡って、evalのコードを見るのが娯楽ですよ。やや、さりげなく JIT なんて語句が見えるぞ。bytecodeを更にマシン語に落とし込んで、加速してるのか。 こういう裏側には、普通のユーザーは関係ないだろうけど。
PDF形式のファイルの用紙サイズを調べる方法は？
PDFファイルを作っているんだけど、用紙サイズってどうなるの？ ぐぐったら、アドビの何たらこうたらばかりが出てくる。もうfirefoxで十分でねぇ。しつこく調べたら
PDF形式のファイルをテキストエディタで開き （注：バイナリファイルですので警告画面等が 表示されることがあります） ， MediaBox という文字列を探します．この属性はページの 大きさ （印字可能範囲） を示しており，A4縦の場合は通常 [0 0 595 842] が設定されて います （単位は1/72インチ = 1ポイント）．
sakae@ub:~/ANA$ grep MediaBox ZZZZ.pdf Binary file ZZZZ.pdf matches
本当に調べてくれているのか？ grepはファイルの頭の部分を見て、諦めている風。-a オプションを付けて、強引にテキストファイルと解釈させればおｋ。
sakae@ub:~/ANA$ grep -a MediaBox ZZZZ.pdf << /Type /Pages /Kids [ 7 0 R ] /Count 1 /MediaBox [0 0 841 595] >>
これに気が付く前は、
sakae@ub:~/ANA$ strings ZZZZ.pdf | grep Media << /Type /Pages /Kids [ 7 0 R ] /Count 1 /MediaBox [0 0 841 595] >>
bainayの海を平気で渡って行く strings なんてのと併用してた。バイナリアンである。
上記は、
pdf('ZZZZ.pdf', paper='a4r' )
として、90度のローテをさせた作成である。