ggplot2で任意の線分

R

baseライブラリ中のlines()に対応するggplot2のレイヤー関数はgeom_line()であると説明される場合が多く、大抵の場合はそれで事足りるのだが、任意の位置に線分を描き込みたいという場合はあまり使い勝手が良くない。
例えば、「統計学：Rを用いた入門書（Michael J. Crawley）」の170ページに載っている次のグラフは、lines()を使ってデータポイントから平均値への線を書いている（詳しくは書籍を参照）。
f:id:Rion778:20151212221826p:plain
同じようなことをggplot2を使ってやろうとして少し調べたところ、結局こうなった。

library(ggplot2)
library(scales)

gaTRUE <- oneway$garden=="A"
gbTRUE <- oneway$garden=="B"

ggplot(oneway, aes(x = 1:length(ozone), y = ozone, color = garden)) +
  geom_point() + 
  geom_hline(y=mean(oneway$ozone[gaTRUE]), color = hue_pal()(2)[1]) +
  geom_hline(y=mean(oneway$ozone[gbTRUE]), color = hue_pal()(2)[2]) +
  geom_segment(aes(x=1:length(ozone), y=ozone,
                   xend=1:length(ozone), 
                   yend=ifelse(gaTRUE,
                               mean(oneway$ozone[gaTRUE]),
                               mean(oneway$ozone[gbTRUE])))) +
  theme_bw() + theme(panel.grid=element_blank())

f:id:Rion778:20151212221843p:plain
線分を引くレイヤー関数としてはgeom_segment()が用意されているので、これを使えば良い。
直接値を指定してgeom_segment(x=1, y=1, xend=2, yend=2)などという使い方も出来る。
ついでに、abline(h=...)に相当するレイヤー関数はgeom_hline()であるが、色の指定をデフォルトのカラーパレットに合わせようとするならばscalesパッケージのhue_pal()をこのように使う。

参考

分布の集中度を調べる

R

ランダムか、集中か

水田の調査圃場において、正方形に区切った調査区画毎にニカメイガの卵塊が次のように見つかったとする。

egg <- c(2,4,1,1,
         4,0,3,3,
         4,3,0,1,
         3,2,1,1,
         2,1,4,3,
         2,2,1,0,
         1,0,5,1,
         2,1,3,2,
         2,2,2,2,
         4,1,1,2,
         3,2,0,1,
         1,2,1,0)

もし、ニカメイガの産卵が先に卵塊があるか否かに影響されず、ランダムに行われるのであれば、その分布は二項分布に従うことが予想される。
また、今回の例のように、任意の個体の特定の調査区への飛来確率pが非常に小さく*1、総個体数Nが大きい場合、二項分布をポアソン分布で近似できる。
よって、得られた分布のポアソン分布への適合度を確認することで分布がランダムか否かを確認できる。
また、ポアソン分布は分散と平均値が等しいという性質を持つことから、得られた分布の分散と平均の比\sigma^2/m（実際にはs^2/\bar{x}）が1に等しいか、大きいか、小さいかを確認することでもポアソン分布かどうかの判断ができる。

ポアソン分布への適合度を確認する

まず、得られた分布から度数分布表を作成する。そして、分布がポアソン分布に従っていた場合の各度数の理論値を求める。ポアソン分布はパラメータとして母平均mのみを持つ。ここでは平均値\bar{x}を母平均の不偏推定量として用いる。

> ## ランダム分布か、集中分布か（ポアソン分布への当てはめ）
> # 度数分布の計算
> (e.table <- table(egg))
egg
 0  1  2  3  4  5 
 6 15 14  7  5  1 
> # ポアソン分布における各度数の比率計算
> (e.table <- data.frame(count = c(e.table, 0, 0),
+ p = dpois(0:7, mean(egg))))
  count           p
1     6 0.156583374
2    15 0.290331673
3    14 0.269161656
4     7 0.166356857
5     5 0.077113335
6     1 0.028596195
7     0 0.008837019
8     0 0.002340758
> # 最後の項はまとめる
> e.table$p[8] <- 1 - sum(e.table$p[1:7])

比率は合計が1になるように、分布に現れなかった項まで計算し、余分な分はまとめておく。

> # 区画数の理論値計算
> (e.table <- cbind(e.table, phi = e.table[,2] * length(egg)))
  count           p        phi
1     6 0.156583374  7.5160020
2    15 0.290331673 13.9359203
3    14 0.269161656 12.9197595
4     7 0.166356857  7.9851291
5     5 0.077113335  3.7014401
6     1 0.028596195  1.3726174
7     0 0.008837019  0.4241769
8     0 0.003019892  0.1449548
> # 理論値5以下の区画はまとめる
> (e.table <- rbind(e.table[1:4,], lapply(e.table[5:8,], sum)))
  count         p       phi
1     6 0.1565834  7.516002
2    15 0.2903317 13.935920
3    14 0.2691617 12.919759
4     7 0.1663569  7.985129
5     6 0.1175664  5.643189

ここで理論値が小さな項目をまとめるのは、観測値は0と1の間をとらないため、後にカイ二乗統計量を計算する際にこれが大きくなりすぎるのを防ぐため、ということのようである。
次にカイ二乗統計量＝観測値と理論値の差の自乗の総和を求める。

> # 理論値との差
> (e.table <- cbind(e.table, d = abs(e.table$count - e.table$phi)))
  count         p       phi         d
1     6 0.1565834  7.516002 1.5160020
2    15 0.2903317 13.935920 1.0640797
3    14 0.2691617 12.919759 1.0802405
4     7 0.1663569  7.985129 0.9851291
5     6 0.1175664  5.643189 0.3568109
> # 差の自乗
> (e.table <- cbind(e.table, d2 = e.table$d^2))
  count         p       phi         d        d2
1     6 0.1565834  7.516002 1.5160020 2.2982620
2    15 0.2903317 13.935920 1.0640797 1.1322655
3    14 0.2691617 12.919759 1.0802405 1.1669196
4     7 0.1663569  7.985129 0.9851291 0.9704794
5     6 0.1175664  5.643189 0.3568109 0.1273140
> # カイ二乗統計量の計算
> (e.table <- cbind(e.table, "d2/phi" = e.table$d2/e.table$phi))
  count         p       phi         d        d2     d2/phi
1     6 0.1565834  7.516002 1.5160020 2.2982620 0.30578251
2    15 0.2903317 13.935920 1.0640797 1.1322655 0.08124799
3    14 0.2691617 12.919759 1.0802405 1.1669196 0.09032054
4     7 0.1663569  7.985129 0.9851291 0.9704794 0.12153584
5     6 0.1175664  5.643189 0.3568109 0.1273140 0.02256065
> (chi <- sum(e.table$"d2/phi"))
[1] 0.6214475

そして得られたカイ二乗統計量を検定にかける。自由度は階級数-2=5となる。

> # 検定（自由度=階級数-2）
> pchisq(chi, 3, lower.tail = FALSE)
[1] 0.8915054

以上より分布がポアソン分布に従わないとは言えない。
ニカメイガの卵塊はポアソン分布に近い分布をするだろう（実際にはプロットしてみたほうが良い）。

分布の集中度を調べる

上述のように分布がポアソン分布に従うのであれば、s^2/\bar{x}を調べることで分布の集中度、一様度が判定出来る。s^2/\bar{x}が

  • 1より小さい = 一様分布
  • 1である = ランダム分布（ポアソン分布）
  • 1より大きい = 集中分布

である。
さらに、s^2/\bar{x}はF分布するので、F検定によって分布が集中分布であるか否かを検定できる。自由度はn_1=n-1, n_2=\inftyを用いる。
s^2/\bar{x}が1より小さい場合は、\bar{x}/s^2について検定を行うことで一様分布であるか否かの検定をする。

> ## バリアンスと平均値の関係を確認する
> var(egg)/mean(egg)
[1] 0.8489123
> qf(0.95, length(egg)-1, Inf)
[1] 1.361726
> pf(var(egg)/mean(egg), length(egg)-1, Inf, lower.tail=FALSE)
[1] 0.7590387
> pf(mean(egg)/var(egg), length(egg)-1, Inf, lower.tail=FALSE)
[1] 0.188295

ここで注意が必要なのは、s^2/\bar{x}は平均値によって値が変わってしまうため、s^2/\bar{x}が大きいほど分布が集中している、とは言えないという点である。
この点を改良した指数として、I_\delta指数（Morisita's indexとも）がある。この指数はs^2/\bar{x}と同様にポアソン分布で1、一様分布で1より小、集中分布で1より大となるが、平均値に影響されず*2、平均値の異なる集団同士の集中度の比較に用いる事ができる。
I_\delta=\frac{\sum x_i(x_i-1)}{N(N-1)}
ここでx_iはi番目の区画における個体数、Nは総個体数である。

> ## I_δ示数の利用
> i.delta <- function(x){
+ length(x) * (sum(x * (x - 1))/(sum(x)*(sum(x)-1)))
+ }
> i.delta(egg)
[1] 0.9193054

参考文献

*1:二項分布の平均値はNpに等しいことから求めると、約2%

*2:一部例外はある

Rで逆推定

R

JMPには逆推定というコマンドがある

次のようなデータを用意して、線形回帰を行うとする。

temp <- rep(c(15, 20, 25), c(5, 5, 5))
day <- c(34, 33, 36, 37, 35,
         30, 28, 29, 25, 28,
         23, 25, 22, 24, 26)

こんなかんじで適当に。

> result <- lm(day~temp)
> plot(day~temp)
> abline(result)
> summary(result)

Call:
lm(formula = day ~ temp)

Residuals:
   Min     1Q Median     3Q    Max 
 -4.00  -1.00   0.00   1.25   2.50 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   51.000      2.307  22.110 1.07e-11 ***
temp          -1.100      0.113  -9.734 2.46e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.787 on 13 degrees of freedom
Multiple R-squared:  0.8794,	Adjusted R-squared:  0.8701 
F-statistic: 94.76 on 1 and 13 DF,  p-value: 2.457e-07

f:id:Rion778:20150805235054p:plain
ここで、統計ソフトのJMPで同じようなことをしたとすると、この後に逆推定というやつができる。
すなわち、上記の例で言えば、dayが特定の値（例えば30）になるようなtempの値を求めるという操作である。このとき、tempの信頼区間も同時に出力される。
ちなみに、上記データに基づいてJMPでtemp=30について信頼水準=0.95で逆推定を行うと、

  • 予測値: 19.09091
  • 下側95%:18.09069
  • 上側95%:19.99694

という推定値が出力される。
どうもこれがRではそんなに簡単にはできないらしい。関数やパッケージの調査不足かもしれないので、良いやり方をご存じの方があれば教えて頂きたい*1。

回帰式を変形させる

y = a + bx
という回帰式を、
y = A + b(x-c)
という風に変形する。このとき、Aは目的とするyの値（上記例で言えば30）であるとすると、x-c=0のときにy=Aとなるので、cが求めるべきxの値に等しいとういことになる。あとは、bとcについて、非線形回帰を用いて求め、cの信頼区間を求める。
（参考）

非線形回帰および信頼区間

非線形回帰はnlsを用いて行う。startの値はどう選ぶのが適切かはよく分からないが、選び方がまずいと動かない時がある。

result2 <- nls(day~30+b*(temp-c), start=c(b=0.1,c=0))

予測値はnlsオブジェクトの中身を見れば書いてある（これは回帰直線上の話だからわざわざ非線形回帰をしなくても分かるが）。

> result2
Nonlinear regression model
  model: day ~ 30 + b * (temp - c)
   data: parent.frame()
    b     c 
-1.10 19.09 
 residual sum-of-squares: 41.5

Number of iterations to convergence: 8 
Achieved convergence tolerance: 1.162e-09

パラメーターの信頼区間は、confintで求める。

> confint(result2)
Waiting for profiling to be done...
       2.5%      97.5%
b -1.344124 -0.8558761
c 18.090557 19.9969627

JMPの結果とほんの少し違うのが気になるところ。

*1:Rではlmオブジェクトに対してpredictという関数が用意されているが、これは例で言えばtempが特定の値のときのdayの推定値と信頼区間を求めるというもので、逆推定はできない。また、chemCalというパッケージの中にinverse.predictといういかにもな名前の関数が入っているが、それではイマイチ上手く行かなかった。使い方が間違っているのかもしれない。

不審なファイルの存在確認等

Windows

事の発端

www.ipa.go.jp
「バッチでの配布でも可能です」とか書くならバッチファイル用意しといてほしい…。

バッチファイル

上からやれといわれたけど手順書の通りにはとてもやってられないので渋々バッチファイルを作成した。
こいつが必要な立場にある人はgithubとかアップローダーとかにアクセスできないと思うので直接書いておく。メモ帳か何かにコピペしてhoge.batとか適当な名前+.batで保存して実行のこと。
ちなみにWindows XP以前ではスタートアップフォルダの位置が異なるのでそのままでは動作しない。スクリプト中のスタートアップフォルダ記述箇所を修正すれば動作する。ただし、schtasksコマンドはProfessionalでないと無いそうだ。まあ必要になること自体がおかしいので修正はしない。

@echo off
setlocal enabledelayedexpansion

echo 1. ファイルの存在有無の確認 ---------------------------------
set RESULT=
for /f "usebackq" %%i in (`dir /a /r /s /b "%TEMP%" "%SystemDrive%\Users\%USERNAME%\AppData\Roaming\Microsoft\Windows\Start Menu\Programs\Startup\" "%SystemDrive%\Users\All Users\Microsoft\Windows\Start Menu\Programs\Startup" ^|findstr /I /R "\\leanp\.exe \\leassap\.exe \\leassaq\.exe \\leassnp\.exe \\mdm\.exe \\nvsvcv\.exe \\nvvscv\.exe \\slwga\.exe \\upsl\.dll \\vmat\.exe \\vmatam\.exe \\vmatap\.exe \\vmater\.exe \\vmmat\.exe \\vmnatam\.exe \\vmwere\.exe \\windump\.exe \\ct\.exe \\yrar\.exe \\csvde\.exe \\GetPassword\.exe \\mimikatz\.exe \\mimikatzx64\.exe \\mimikatz1\.exe \\gp\.exe \\Gp64\.exe \\ps\.txt \\msver\.exe \\ss\.exe \\mailfinal\.exe \\mail_noArgv_final\.exe \\result\.log \\14068\.rar \\ms14-068\.exe \\kptl\.doc \\kenpo\.doc"`) do (set RESULT=!RESULT!^

%%i
)

echo;
if not "!RESULT!"=="" (
  echo 下記のファイルが発見されました。
  echo !RESULT!
) else (
  echo 不審なファイルは発見されませんでした。
)
echo;

echo 2-(1) 自動起動設定に不審なファイルが登録されていないかの確認 --
set RESULT=
for /f "usebackq" %%i in (`schtasks /query /v ^|findstr /I /R "\\leanp\.exe \\leassap\.exe \\leassaq\.exe \\leassnp\.exe \\mdm\.exe \\nvsvcv\.exe \\nvvscv\.exe \\slwga\.exe \\upsl\.dll \\vmat\.exe \\vmatam\.exe \\vmatap\.exe \\vmater\.exe \\vmmat\.exe \\vmnatam\.exe \\vmwere\.exe \\windump\.exe \\ct\.exe \\yrar\.exe \\csvde\.exe \\GetPassword\.exe \\mimikatz\.exe \\mimikatzx64\.exe \\mimikatz1\.exe \\gp\.exe \\Gp64\.exe \\ps\.txt \\msver\.exe \\ss\.exe \\mailfinal\.exe \\mail_noArgv_final\.exe \\result\.log \\14068\.rar \\ms14-068\.exe \\kptl\.doc \\kenpo\.doc"`) do (set RESULT=!RESULT!^

%%i
)

@echo off

echo;
if not "!RESULT!"=="" (
  echo 下記のファイルが発見されました。
  echo !RESULT!
) else (
  echo 不審なファイルは発見されませんでした。
)
echo;

echo 2-(2) レジストリキーに不審なファイルが登録されていないかの確認
set RESULT=
for /f "usebackq" %%i in (`reg query "HKLM\Software\Microsoft\Windows\CurrentVersion\Run" /s ^|findstr /I /R "\\leanp\.exe \\leassap\.exe \\leassaq\.exe \\leassnp\.exe \\mdm\.exe \\nvsvcv\.exe \\nvvscv\.exe \\slwga\.exe \\upsl\.dll \\vmat\.exe \\vmatam\.exe \\vmatap\.exe \\vmater\.exe \\vmmat\.exe \\vmnatam\.exe \\vmwere\.exe \\windump\.exe \\ct\.exe \\yrar\.exe \\csvde\.exe \\GetPassword\.exe \\mimikatz\.exe \\mimikatzx64\.exe \\mimikatz1\.exe \\gp\.exe \\Gp64\.exe \\ps\.txt \\msver\.exe \\ss\.exe \\mailfinal\.exe \\mail_noArgv_final\.exe \\result\.log \\14068\.rar \\ms14-068\.exe \\kptl\.doc \\kenpo\.doc"`) do (set RESULT=!RESULT!^

%%i
)

@echo off

echo;
if not "!RESULT!"=="" (
  echo 下記のファイルが発見されました。
  echo !RESULT!
) else (
  echo 不審なファイルは発見されませんでした。
)
echo;

endlocal

echo;
pause

調べたこと

簡単に出来るかと思ったら結構面倒くさくて色々と調べる羽目になった。

手順書見るとTEMPフォルダ2回確認してるし、拡張子が違う場合があると書いてあるのに記載のコマンドでは.exeしか見てないし色々とアレ。

ggplot2で縦に並べたグラフの横幅を揃える

R

（※2015/06/03 20:27追記あり）
下記の記事に基づく発言に関連してどうも某所で勘違いが発生しているようなので。www.exegetic.biz

使用データ

下記のように生成したものを用いる。

x <- 0:100
x <- data.frame(x = x, y1 = sin(x * pi / 10), y2 = x^2)

なぜ横幅がずれるのか

まず、これをbaseで普通にプロットする。

layout(matrix(1:2, 2))
plot(y1 ~ x, type = "l", data = x)
plot(y2 ~ x, type = "h", data = x)
par(old)

f:id:Rion778:20150602235943p:plain
グラフの横幅はずれない。
baseを用いたプロットでは、あえてずらすような細工をしない限り、「グラフの横幅のずれ」は発生しない。
よって、縦軸を揃えるためにマージン等を設定する必要は無い（そのままでは余白が大きくなりすぎるという問題はあるが）。
ちなみに、y軸の目盛りラベルが話にならない*1ので、las = 1を設定して再確認してみるもやはりずれない。

layout(matrix(1:2, 2))
plot(y1 ~ x, type = "l", data = x, las = 1)
plot(y2 ~ x, type = "h", data = x, las = 1)
par(old)

f:id:Rion778:20150603000700p:plain
ずれないが、下のグラフのy軸ラベルと目盛りラベルが重なってしまっている。
といったところでピンとくるかもしれないが、gridの作図では*2、baseの作図と以下の点が異なる。

  • y軸の目盛ラベルは水平方向がデフォルト（las = 1を指定した場合と同様）。
  • y軸の目盛ラベルの幅に応じて余白が自動調整される。

これは通常のプロットではほとんど問題にならないが、y軸ラベルの長さが異なるグラフを縦に並べると、調整なしでは横幅が必ずずれるという弊害になって現れる。

library(ggplot2)
library(gridExtra) # grid.arrange()関数のために必要
p1 <- ggplot(x, aes(x = x)) + geom_line(aes(y = y1)) + theme_classic()
p2 <- ggplot(x, aes(x = x)) + geom_bar(aes(y = y2), stat = "identity") + theme_classic()
grid.arrange(p1, p2)

f:id:Rion778:20150603001835p:plain
gridの作図ではlayoutを用いたレイアウトはできないので、grid.arrangeなどの関数を用いる必要がある*3。
なお、最初に提示したリンク先の例では上下のグラフの高さも設定しているので、これとは出力が異なっている。

解決策

gridの作図では、作図オブジェクト（graphical object; grob）というオブジェクトが生成され、これを通じて描画の制御ができる。
よって、横幅を揃えるためには、まずgrobから横幅や縦幅といったレイアウトの情報を取り出す*4。

g1 <- ggplot_gtable(ggplot_build(p1)) 
g2 <- ggplot_gtable(ggplot_build(p2))

ggplot_build関数は描画に必要な全ての情報をplotオブジェクトから生成し、ggplot_gtable関数はそこから横幅や縦幅といったレイアウトに必要な情報を含むgtableオブジェクトを生成する*5。
gtableオブジェクトから横幅の情報を取り出し、これをp1とp2のどちらかに合わせる。合わせる基準としては、幅の広い方を選択することにする（文字の重なりを防ぐため）。
gridパッケージには単位を含むオブジェクトの要素ごとに大小を比較するunit.pmax、unit.pmin関数が含まれているので、これを用いて横幅の値の大きい方を取り出す。

maxWidth <- unit.pmax(g1$widths, g2$widths)

この情報を用いてグラフのパラメータを上書きすることで、横幅を揃えることができる。

g1$widths <- maxWidth
g2$widths <- maxWidth
grid.arrange(g1, g2)

f:id:Rion778:20150603003838p:plain
要するに、gridに起因する横幅のずれの調整にggplot2に含まれる関数が使えるという話。

追記


ということを教えていただいた。
size=に"max"や"min"を指定するとエラーが出るが、gtableをこれhttps://github.com/baptiste/gtableに差し替えれば動作する様子。first/last指定でも大抵は事足りるとは思うが。

*1:そもそもどうしてこの方向がデフォルトなのか今ひとつ納得出来ない。

*2:つまるところggplot2が元凶ではない。

*3:grid.arrangeはarrangeGrob関数のラッパーで、arrangeGrob関数はかなり色々なカスタマイズが便利にできる関数の様だが話がずれるのであまり調べていない。

*4:ggplot2はこの2つの関数を含むので、gridに起因する横幅問題の解決が容易になる、というのが元記事の要点である。

*5:と思うがこれもあまり良く調べていない。

RStudioの折りたたみ機能

R

今更ながらコードの折りたたみ機能が付いている事に気付いた。
詳細はCode Folding and Sections – RStudio Supportを参照。

折りたたみ設定されるもの

  • ブレース（波括弧{}）で括られた領域。
  • R SweaveやR Markdownドキュメントにおけるコードチャンク。
  • R Mardcownドキュメント中におけるセクション。
  • コードセクション（後述）

上記4つは自動的に折りたたみ設定がされる。
例えば、関数や条件定義領域などのブレースで囲まれた部分があると、行数の横に三角が現れる。
f:id:Rion778:20150531170144p:plain
これをクリックすることで折りたたみができる。
f:id:Rion778:20150531170225p:plain
展開する場合は再び三角をクリックするか、折りたたまれている事を示している青いアイコンをクリックする。
（キーボードショートカットで操作する方法は用意されていないか、うまく動かないようだ。Ctrl + Command + L（Mac）で閉じられるときもあるし、閉じられないときもある。このコマンドは任意領域の折りたたみの際には上手く動作する。）

任意の折りたたみ領域を設定する

折りたたみ設定がされていない領域も折りたたむことができる。
範囲選択後にEdit > Folding > CollapseまたはCtrl + Command + L（mac）で任意領域の折りたたみができる。
f:id:Rion778:20150531171014p:plain
f:id:Rion778:20150531171158p:plain

コードセクション

Rスクリプト中に、コメント記号（#）で開始し、4つ以上の-または=または#で終了する行があると、セクションの開始とみなされる。
f:id:Rion778:20150531171605p:plain
どのセクション記号も機能は同様であり、セクション開始から次のセクション開始までが折りたたみの対象となる（＝サブセクションを作ることはできない）。
セクションは折りたたみの対象となるだけでなく、関数定義と同様にジャンプ機能の対象となり、その部分へ容易に移動できるようになる。
f:id:Rion778:20150531171944p:plain

キーボードショートカット

Macの場合のショートカットが若干押しにくいが…

機能 Windows Mac
折りたたみ Alt + L Ctrl + Command + L
展開 Alt + Shift + L Ctrl + Command + Shift + L
全て折りたたみ Alt + O Ctrl + Command + O
全て展開 Alt + Shift + O Ctrl + Command + Shift + O

Edit > Foldingから選択することも出来る。
折りたたみ・展開は基本的には領域選択と合わせて使う。領域を選択していない場合は動くこともあるし動かないこともあってよく分からない。
「全て折りたたみ」を実行した場合、最も「外側」で折りたたまれる。もしセクションが設定されていればセクションが最も外側となる。
なお、折りたたみの状態はファイルを閉じない限り保存される（たとえRStudioを終了しても）が、ファイルを閉じて再度開いた際には全ての折りたたみが展開された状態となる。

S3

R

RStudioではじめるRプログラミング入門を読んだ。RStudioについて詳しい情報があるかと思っていたのだが、これについてはそこそこの記述だった。
ただ、全体的に見て、Rとプログラミングをゼロから学ぼうと思うのであれば良い情報量だと感じた。リゲス本なんかに比べるとずっと表現が平易で分かりやすい。タイトルに「入門」と付いたオライリーの本が入門書であったことなど殆ど無いので逆に面食らった。ゆえに、本当の初心者で無ければあえて読む必要も無い本とも言えるのだが、環境やクラス、コードのベクトル化など、なかなか初心者向けの丁寧な解説が無い部分の記述もあり、個人的には参考になった。また、「他人にRを教える」場合の参考書としても良いのではないかと思う。

S3に関してよくまとまっていたので覚書。なお、本書の説明のほうがずっと長いが丁寧。

S3とは

Rに備わっているクラスシステムの一つで、ジェネリック関数、メソッド、クラスに基いている。

やること

  1. オブジェクトのクラス名を考える
  2. クラスのインスタンスのclass属性にクラス名を設定する
  3. 必要となりそうなジェネリック関数のクラスメソッドを書く

ジェネリック関数

Rには、与えられたオブジェクトの種類によって異なる振る舞いをする関数が備わっている。print()はその例。

> x <- 100000
> print(x)
[1] 1e+05
> attr(x, "class") <- "POSIXct"
> print(x)
[1] "1970-01-02 12:46:40 JST"

同じprint(x)でも"class"属性が変わることで関数の動作が変わっている。このように、ジェネリック関数はクラスに応じて振る舞いを変える。

メソッド

ジェネリック関数は呼び出された際にUseMethod()という特別な関数を呼び出す。

> prin
function (x, ...) 
UseMethod("print")
<bytecode: 0x10432ebc8>
<environment: namespace:base>

UseMethod()は関数の第1引数に与えられたオブジェクトのクラスを調べ、そのクラスを処理するための関数に全ての引数を渡す。
具体的には、print.POSIXct()の様に関数名とクラス名をドットで区切った名称の関数が呼ばれる。

> print.POSIXct
function (x, ...) 
{
    max.print <- getOption("max.print", 9999L)
    if (max.print < length(x)) {
        print(format(x[seq_len(max.print)], usetz = TRUE), ...)
        cat(" [ reached getOption(\"max.print\") -- omitted", 
            length(x) - max.print, "entries ]\n")
    }
    else print(format(x, usetz = TRUE), ...)
    invisible(x)
}
<bytecode: 0x10c0b90e8>
<environment: namespace:base>

これらの関数は（print()の場合）、printのメソッドと呼ばれる。どのようなメソッドが存在するのかについては、methods(print)のようにして調べることができる。

> methods(print)
  [1] print.abbrev*                                
  [2] print.acf*                                   
  [3] print.anova*                                 
  [4] print.Anova* 
...

自作のクラスを追加する

オブジェクトにクラス属性を追加する方法はいくつかある。attr()を使う方法や、class()を使う方法、structure()を使う方法などである。

x <- 1

class(x) <- "myclass"
attr(x, "class") <- "myclass"
structure(x, class = "myclass")

structure()を使うやり方は関数の返り値を構成する場合などに役立つだろう。

> roll <- function(){
+   # サイコロを2つ振って合計値を求める
+   die <- 1:6
+   dice <- sample(die, size = 2, replace = TRUE)
+   structure(sum(dice), dice = dice, class = "dice")
+ }
> 
> roll()
[1] 10
attr(,"dice")
[1] 6 4
attr(,"class")
[1] "dice"

class属性はメソッドがなければ単なる属性扱いであり、デフォルトのprint()では他の属性と同様に表示される*1。

自作のメソッドを作成する

自作のメソッドを作成、追加する方法はとても簡単で、処理方法を定めた関数を作成し、他のメソッドと同様に、関数名とクラス名をドットで区切った関数名を設定すれば良い。

> print.dice <- function(x, ...){
+   dice <- attr(x, "dice")
+   dice <- paste("dice1:", dice[1], " dice2:", dice[2])
+   string <- paste(x, dice, sep = "\n")
+   cat(string)
+ }
> 
> roll()
11
dice1: 6  dice2: 5

*1:printはコンソール上では暗黙に使われるので上記例では明示的に使用していない。

S-Insertでペーストできない

Vimperator Mac

.vimperatorrc中でペースト操作をしたいところに今までS-Insertと書いていたのだけど、気付いたら動作しなくなっていた。
今まではVimperator側でペースト操作に上書きしていたが（vimperatorでのコマンドラインへの貼り付け方法が分からなかった, S-Fマガジン 2008年 04月号 [雑誌], 近所のホームセンターで前飼っていた犬に、そっくりな犬がいた - ぽっぺん日記@karashi.org(2008-04-12)）、Vimperator 3.9から

We don't override ＜S-Insert＞ in textboxes anymore,

とかなってるのでこれが原因だとは思う。
S-InsertでペーストできないのがOS Xの正しい動作なのか、Karabinerとかで色々いじってるのが原因なのかはよく分からない。
とりあえず、S-InsertはM-vに置き換えた。

■

Agriculture Diary

キュウリに対して、おそらく同時期に放飼したであろうスワルスキーカブリダニが容易に見つかるほ場と30分探しても見つからないほ場があった。
スワルスキーカブリダニはナスやキュウリの害虫を食べるダニで、生き物だが農薬として登録されている。生き物だろうが天然物だろうが、農作物に対して病害虫や雑草の駆除目的で使って良いのは農薬だけであって、農薬以外のもの=安全性の担保されていないものを防除目的で使用してはいけない。そして、登録通りに放飼したとすると1m^2あたり25〜50頭導入される。
施設のキュウリというのは1m^2にだいたい1本植わっていて、主枝摘心のつる下ろしという栽培方法だと1株から4〜5本のつるが伸びていて、1本のつるには少なくとも14〜15枚、摘葉をやめていれば20枚かそれ以上の葉がついている。そうすると1m^2に50枚とか100枚とかいう葉があって、そこに25〜50頭なので、よく分散しているとして2〜3枚に1頭は見つかって欲しい、という気になる。
これが見つからないとすると、考えられる理由は、

  1. スワルスキーカブリダニを殺してしまう農薬を使ってしまった。
  2. 増殖するために充分な温度が確保できていない。
  3. 増殖するために充分な湿度が確保できていない。
  4. 気が早い。

1から3は割とありがちで、1ヶ月たっても居ないとするといずれかを疑ったほうが良い。
実は4が特にありがちで、スワルスキーカブリダニというのは結構バカにならない値段がするもので、放飼するのも結構時間がかかるものだから、入れて見えないと不安になる。ただとても小さいもので、また別に葉っぱに均等に分散してくれるわけでもないので、最初の数週間どう頑張っても見つからないということは良くある。だから見えなければ（害虫が増えない限り）しばらく待つのが正解だが、もし自分が金払って労力かけて放飼したものだとしたら、やはり気が気でないだろうと思うところ。