ggplot2で凡例のラベルと項目名を操作する

ラベルと項目名を操作する

以下の説明ではggplot2を含めて3つのパッケージを使う。

library(dplyr)
library(tidyr)
library(ggplot2)

データは次のように準備した。

# テストデータの準備
testdata <- data.frame(x = seq(1, 10, 0.1)) %>%
  mutate(sin = sin(x), cos = cos(x))

散布図形式のプロットに適した形にデータフレームを整形する。

# 散布図の骨格を決める
p <- testdata %>%
  gather(ftype, val, -x) %>% # x以外の変数をvalにまとめ、列名をftypeにグループ変数として格納する
  ggplot(aes(x = x, y = val))

geom_line()を使ってグラフを描く。このとき、colorやltyなどの審美的属性にグループ変数をマッピングすることで変数をグループ別にプロットできる。そして、凡例は勝手に追加される*1。

# colorをftypeにマッピングする
p + geom_line(aes(color = ftype))

凡例のラベルを操作する簡単な方法はlabs()関数を使うこと*2で、例えば

p + geom_line(aes(color = ftype)) + labs(color = "関数")

のようにしてやればよい。
しかし、項目名（ここでのsinとcos）まで変更しようとすると、labs関数では設定ができない。項目名を設定するには、スケールを調整する関数を使用する必要がある。カラースケールの場合はscale_color_◯◯()関数にlabels引数を設定する。◯◯の部分にはhueやgreyなどスケール名が入るが、デフォルトの配色ではhueである。

# colorをftypeにマッピングし、凡例のラベルを変更する
p + geom_line(aes(color = ftype)) +
  scale_color_hue(name = "関数", labels = c(sin = "正弦", cos ="余弦") ) +
  theme_grey(base_family = "HiraKakuProN-W3") # OS Xで日本語をプロットする際にはbase_familyの指定が必要

複数の審美的属性が1つのグループ変数に割り当てられている場合

1つのグループ変数に複数の審美的属性を同時にマッピングすることもできる。

# colorとltyの2属性にftypeをマッピングしたプロットを作成
p2 <- p + geom_line(aes(color = ftype, lty = ftype))
# そのままプロット
p2

このとき、ラベルの片方だけを変更してしまうと、凡例が2つに分裂してしまうので注意が必要である。

theme_set(theme_grey(base_family = "HiraKakuProN-W3")) # themeのデフォルト設定そのものを変更
# 凡例のラベルをcolorの方だけいじった場合
p2 + scale_color_hue(name = "関数", labels = c(sin = "正弦", cos ="余弦"))

これを防ぐには、2つの凡例に全く同じラベル名と項目名を設定する。なお、線種のスケールはscale_linetype()で調整できる。

# 凡例のラベルを両方いじった場合
p2 + scale_color_hue(name = "関数", labels = c(sin = "正弦", cos = "余弦")) +
  scale_linetype(name = "関数", labels = c(sin = "正弦", cos = "余弦"))

labs()を使ってラベルだけを変更する場合も同様で、

p2 + labs(color = "関数", lty = "関数")

のようにしてやる必要がある。

*1:baseパッケージでは凡例は基本的にはlegend()関数などで明示的に追記する必要がある。

*2:labs()関数は他にもtitleやx軸/y軸ラベルの設定もできる。

facetのstripのtextを変更する

ggplot2でfacet_gridやfacet_wrapを使ってグラフを分割描画した際に、各グラフについているラベル（strip）のテキストは、何もしなければ分割に用いたグループ名が自動的に入る。
例：

ggplot(iris, aes(x = Petal.Length, y = Petal.Width)) +
  geom_point() + facet_grid(.~Species)

このラベル部分はlabeller functionという関数で制御されていて、標準でいくつかの関数が用意されている（help(labellers)を参照のこと）。
これをマニュアルで設定するためには、これまではlabeller functionを自分で書く必要があった（cf. r - ggplot: How to change facet labels? - Stack Overflow)。定義したlabeller functionをfacet_gridの引数labellerに与えることで目的が達成される。
例：

# これまでのlabeller functionの作り方
labeli <- function(variable, value){
  names_li <- list("setosa" = "せとさ",
                   "versicolor" = "ばーじからー",
                   "virginica" = "ばーじにか")
  return(names_li[value])
}

theme_set(theme_gray(base_family="HiraKakuProN-W3"))
ggplot(iris, aes(x = Petal.Length, y = Petal.Width)) + 
  geom_point() +
  facet_grid(.~Species, labeller = labeli)

ggplot2のバージョンアップにともなって、labeller functionをもう少し分かりやすく書ける関数、as_labellerが追加された。

labeli <- as_labeller(c(`setosa` = "せとさ",
                         `versicolor` = "ばーじからー",
                         `virginica` = "ばーじにか"))

出来上がったlabeller functionの使い方に変更は無い。

飽和水蒸気圧の求め方

何年か前に飽和水蒸気圧の求め方に関して、Goff-Gratchの式やTetensの式やMurrayの式について少し触れたことがあった*1のだが、最近Wikipediaを見ていたらTetensの式だと思っていたものが

Tetens(1930)のパラメータ値によるAugust他の式

飽和水蒸気量 - Wikipedia

という表記になっていた。この部分はTetensの式となっていたものが割と最近更新された様子である。

Tetensの式はTetensの式ではない

Tetensの式やMurrayの式と呼ばれているこれらの式は、Magnusの式（後述するがこの呼称も適切ではない）と呼ばれている次式、

のパラメータA, B, Cを定めたものに相当し、基本的には同様の式である。Tetensのパラメータがそうであるように、自然対数部分を常用対数として使うパラメータもある。
Tetens(1930)*2によるパラメータを用いたものがTetensの式、Murray(1966)*3のパラメータを用いたものがMurrayの式として呼ばれていたものに該当する。Tetensのパラメータが特に有名となっているのは、少ない桁数でそこそこ正確な近似を与えたという点で優れていたから、ということの様である。

August-Roche-Magnusの式

英語版Wikipediaのクラウジウス-クラペイロンの関係式のページを見ると、飽和水蒸気圧の求め方が応用例としてあげられており、次のような記述がされている。

Fortunately, the August-Roche-Magnus formula provides a very good approximation, using pressure in hPa and temperature in Celsius:

[7][8]

(This is also sometimes called the Magnus or Magnus-Tetens approximation, though this attribution is historically inaccurate.[9])

Clausius–Clapeyron relation - Wikipedia, the free encyclopedia

出典としてあげられているLawrence(2005)*4にもあるように、この式についての初期の仕事はAugust(1828)によってなされており、August-Rocheの式またはAugust-Roche-Magnusの式と呼ぶのがふさわしいようだ。このあたりに基づいて日本語版Wikipediaの記述も変更されたのだろう。

非August-Roche-Magnus式

August-Roche-Magnusの式の形をしていない飽和水蒸気圧を求める式として、Goff-Gratch(1946)の式や、Wexler(1976)の式、Sonntag(1990)の式などがあり、これらはより正確な飽和水蒸気圧を与えるとされている。
計算そのものはコンピュータを用いればすぐに終わるものだが、いずれも式が長く係数の桁数が多いので、そこそこの精度があれば良いという用途には不便である。

結局どうやって求めるの

英語版のクラウジウス-クラペイロンの関係式のページから引用した式の係数は、日本語の教科書によく載っている値や、日本語版Wikipediaに載っている値と若干違う。
この値は、Alduchov and Eskridge(1996)*5によりあたえられたもので、上に述べた3つの非Augst-Roche-Magnus式との相対誤差が、いずれに対しても0.4%以下となる（水上で-40℃〜50℃、氷上で-80℃〜0℃（氷上のものは係数が異なる）の範囲）ように選んだものである。有効な温度範囲も明示されているので、基本的にはこの式を使っておけば良いだろう。
Rの関数として定義する場合の例を示す。

AE1995 <- function(t, s.cool=FALSE){ # 入力：℃　出力：hPa
  ifelse(t > 0 | s.cool,
         6.1094 * exp(17.625 * t / (243.04 + t)), # 水上or過冷却水(s.cool=TRUE)
         6.1121 * exp(22.587 * t / (273.86 + t))  # 氷上
  )
}

Goff-Gratchの式との比較

Goff-Gratchの式をコントロールとして、AE1995を含むいくつかの式の相対誤差をプロットすると、以下のようになる。

日常的な気温の範囲ではどれも大差なく、むしろAE1995は氷点以上のあてはまりがやや悪いようにも見えるが、AE1995はWexler(1976)やSonntag(1990)といった新しい式も含めた評価をしている点を注意しておく必要がある。面倒なので今回はやらないが、これらの式との比較をしてみるのも面白いだろう。
Okadaの式の当てはまりは特に良い様に見えるが、この式はGoff-Gratch式を最小二乗近似したもので、August-Roche-Magnusの式の一種ではない。そのため、元の式よりは簡単であるがやや入力が多い（スクリプトは後に示した）。
氷点以下の飽和水蒸気圧に関しては、AE1995がかなり良い値を示していることが分かるだろう。
なお、Tetensのパラメータは氷点以下のものがよく分からなかったので組み込んでおらず、氷点以下の値は適切ではない。後述のスクリプトを見てもらえばわかるが、Tetensのパラメータはとても簡単で、しかも常用対数を使っている。対象とする温度範囲に氷点以下が含まれないのであればこのパラメータを用いても問題ないだろう。この桁数でこの当てはまりは驚異的であり、Tetensの式と呼ばれてきただけの事はある。

Rのスクリプト

グラフの描画に用いたスクリプトを記しておく。
ライブラリは以下のものを使用した。

library(tidyr)
library(ggplot2)

飽和水蒸気圧を求める式を定義する。AE1995は前述の通り。

GofGra <- function(t, s.cool=FALSE){
  # Goff-Gratchの式
  ifelse((t > 0) | s.cool, 
    10^ # 水または過冷却水(s.cool=TRUE)の場合（-50〜100度）
      (10.79574 * (1 - 273.16/(t + 273.15)) -
         5.02800 * log((t + 273.15)/273.16, 10) +
         1.50475 * 10^{-4} * {1-10^(-8.2969 * ((t + 273.15)/273.16 - 1))} +
         0.42873 * 10^{-3} * {10^(4.76955*(1 - 273.16/(t + 273.15))) - 1} +
         0.78614), 
    10^ # 氷上の場合（-100度まで）
      (-9.09685 * (273.16/(t + 273.15) - 1) -
         3.56654 * (log(273.16/(t + 273.15), 10)) + 
         0.87682 * (1 - (t + 273.15)/273.16) + 
         0.78614)
    )
}

Okada <- function(t, s.cool=FALSE){
  # 岡田の式（Goff-Gratch式を最小二乗近似したもの）
  ifelse((t > 0) | s.cool,
    exp(1.809378 + # 水または過冷却水の場合（-30〜50度）
          0.07266115 * t +
          (-3.003879 * 10^-4) * t^2 +
          (1.181765 * 10^-6) * t^3 +
          (-3.863083 * 10^-9) * t^4),
    exp(1.809378 + # 氷上の場合（-30度まで）
          0.08238957 * t +
          (-2.990908 * 10^-4) * t^2 +
          (1.362765 * 10^-6) * t^3)
  )
}

# August式のパラメータ違い
Murray <- function(t){
  ifelse(t > 0,
         6.1078 * exp(17.2693882 * t /(t + 237.3)),
         6.1078 * exp(21.8745584 * t /(t + 265.5))
         )
  }
Tetens <- function(t){ 6.11 * 10^(7.5*t/(t + 237.3)) }

プロット。tidyrパッケージを使うとプロットに適した形式への変換が簡単だった。

# 対象とする温度範囲
ts <- seq(-100, 100, 0.1)
# 相対誤差を計算しデータフレームを作成
es.diff <- data.frame(air.temp = ts,
                      Okada =  (Okada(ts) - GofGra(ts)) / GofGra(ts) * 100,
                      Murray = (Murray(ts) - GofGra(ts)) / GofGra(ts) * 100,
                      Tetens = (Tetens(ts) - GofGra(ts)) / GofGra(ts) * 100,
                      AE1995 = (AE1995(ts) - GofGra(ts)) / GofGra(ts) * 100)
# longフォーマットへの変換(tidyr使用)
es.diff.l <- gather(es.diff, "method", "diff", -air.temp)
# plot
theme_set(theme_gray(base_family="HiraKakuProN-W3")) # 日本語プロット用設定
ggplot(es.diff.l, aes(x = air.temp, y = diff)) +
  geom_line(aes(color = method)) + xlim(-100, 100) + ylim(-1, 1) +
  labs(title = "Goff-Gratchの式の値に対する相対誤差",
       y = "相対誤差（％）", x = "気温（℃）") 

*1:湿度諸量の計算法(1) - 相対湿度 - - もうカツ丼でいいよな

*2:O. Tetens, "Uber einige meteorologische Begriffe", Z. Geophys., 6, 297–309, 1930

*3:F.W. Murray, "On the Computation of Saturation Vapor Pressure", J. Appl. Meteor., 6, 203-204, 1966

*4:M. G. Lawrence, "The relationship between relative humidity and the dew point temperature in moist air: A simple conversion and applications", Bull. Am. Meteorol. Soc., 86, 225-233, 2005

*5:Oleg A. Alduchov and Robert E. Eskridge, "Improved Magnus Form Approximation of Saturation Vapor Pressure", J. Appl. Meteor., 35, 601–609, 1996

ThinkPad X61 SSD化＋その他いろいろ

2万ほど投資したらまだまだ使えそうな感じになった。

用意したもの

• SSD（A-DATA ASP900S3-256GM-C-7MM ADATA 2.5"SSD 256GB SATA6G A-DATA ASP900S3-256GM-C-7MM）
• 2.5インチHDケース（【日本正規代理店】 ORICO 2.5インチ HDD/SSD 外付け ドライブ ケース SATA対応 USB3.0対応 リムーバブルケース 工具不要 簡単 着脱 人肌のような手触り 2599US3 ブラック）
• メモリ4GB分（Komputerbay 日本進出記念 メモリ 2枚組 DDR2 800MHz PC2-6400 2GBX2 DUAL 200pin SODIMM ノート パソコン用 増設メモリ 2GB デュアル
• Windows7インストールディスク（昔買ったやつ）
• USBフラッシュメモリ（部屋に転がってたやつ）

HDクローン作成

A-DATAのSSDを購入するとAcronis True image HDがダウンロードできるので、これを利用してHDのクローンを作成した。手順としては、

1. Acronis True image HDの起動ディスクを作成
2. シャットダウン（必要であればその後CD/DVDドライブからブートできるようにBIOS設定）
3. SSDをケースに入れてUSB接続
4. 外付けのCDドライブに入れてブート
5. Clone ModeはAutomatic、SouceにHD、DestinationにSSDを選択してスタート

これでSSDにデータがコピーされるので、終了したらHDとSSDを入れ替える。ついでにこのタイミングでメモリも入れ替えた。

Windows7 64bitのインストール

Windows10を入れる前に64bit化をした。Windows7のパッケージ版には32bit版と64bit版両方のインストールディスクが入っているが、このうち64bit版のインストールディスクを入れて新規インストールをした。新規インストールなので、一応データは退避されるが設定などは初期状態に戻る。
ここでWindows Updateを開始したところ全く終わらないので中止してWindows10のインストールを行った。

Windows10のインストール

GetWindows10から無償アップグレードというのが通常の手順だと思うが、これだと途中にWindows Updateが入ってくるのでそこで止まってしまって進まない。

1. https://www.microsoft.com/ja-jp/software-download/windows10の下の方にある「ツールを今すぐダウンロード」からダウンロードできるインストールメディア作成ツールをダウンロード。
2. ツールを使ってUSBフラッシュメモリをインストールメディアにする
3. USBフラッシュメモリからWindows10をインストールする（USBフラッシュメモリからブートするわけではなく、Windows7上で実行する）

このとき、Windows Updateをするかどうか訪ねてくると思うが、これを後回しにすることでインストールを完了できる。

Windows10インストール後の設定

• Metroアプリでのスクロール有効化：Synaptics ThinkPad UltraNav ドライバー (Windows 10 64bit) - ThinkPad E450, E455, E555 - Lenovo Support (JP)からドライバをダウンロードし、デバイスマネージャーから手動でPS/2 Trackpointのドライバのアップデートを行う（参考：

いまさらハック: Win10にアップグレードしたX220のTrackpointセンタースクロールをMetroアプリでも動かす方法）。

• 無線LANが電源復帰時に開始されない問題の修正：デバイスマネージャーからIntel Wireless WiFi Link 1965AGNのプロパティを開き、電源の管理タブにある「電力の節約のために、コンピューターでこのデバイスの電源をオフにできるようにする」のチェックを外す（参考：外部記憶: ThinkPad X61のWindows 10で、無線LANが休止状態から復帰時、接続できない問題）。

Ubuntu再インストール

ここまででWindowsは使えるようになったが、別パーティションに入れていたUbuntuが起動しなくなっていたのでインストールしなおした。

1. Ubuntu Desktop 日本語 Remixのダウンロード | Ubuntu Japanese TeamからUbuntu 15.10のisoファイルを入手。
2. UNetbootin（UNetbootin, Universal Netboot Installer 日本語情報トップページ - OSDN）と上記のisoファイルを利用してUSBフラッシュメモリをインストールメディア化する。
3. インストール（BIOS設定からUSB HDDの起動順位を上げる必要がある）

Windows標準のブートローダーを使用するようにする

ここまで行き当たりばったりでやった結果、起動するとgrubが立ち上がり、grubからWindows10を選択するとWindowsのブートローダーが立ち上がるという面倒な状態になってしまったのでこれを修正する。

• WindowsのブートローダーのOS選択画面にある「規定の起動設定や…」からコマンドプロンプトを起動させる
• 下記コマンドを入力してMBRを修復する

bootrec /fixboot 
bootrec /fixmbr
exit

ここまででgrubは起動しなくなる。次に、Windowsのブートローダーの項目を修正する。

• EasyBCDをインストール（EasyBCD 2.3 - Windows Vista 以降のOS のブートマネージャーを編集するツール）
• EasyBCDを起動し、[Edit Boot Menu]からWindows10以外の項目を一旦消す。
• [Add New Entry]の[Linux/BSD]タブにてTypeをGRUB 2に、DriveをUbuntuがインストールされているパーティションに設定してNameをUbuntuとか適当に設定してAdd New Entryする。

これで起動後はWindows標準のブートローダーが起動するようになる。そこからUbuntuを選択するとgrubが起動する。

ggplot2で任意の線分

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

同じようなことを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())

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

参考

• ggplot2 Quick Reference: geom_segment | Software and Programmer Efficiency Research Group
• r - get list of colors used in a ggplot2 plot? - Stack Overflow

分布の集中度を調べる

ランダムか、集中か

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

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が大きい場合、二項分布をポアソン分布で近似できる。
よって、得られた分布のポアソン分布への適合度を確認することで分布がランダムか否かを確認できる。
また、ポアソン分布は分散と平均値が等しいという性質を持つことから、得られた分布の分散と平均の比（実際には）が1に等しいか、大きいか、小さいかを確認することでもポアソン分布かどうかの判断ができる。

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

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

> ## ランダム分布か、集中分布か（ポアソン分布への当てはめ）
> # 度数分布の計算
> (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

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

分布の集中度を調べる

上述のように分布がポアソン分布に従うのであれば、を調べることで分布の集中度、一様度が判定出来る。が

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

である。
さらに、はF分布するので、F検定によって分布が集中分布であるか否かを検定できる。自由度は, を用いる。
が1より小さい場合は、について検定を行うことで一様分布であるか否かの検定をする。

> ## バリアンスと平均値の関係を確認する
> 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

ここで注意が必要なのは、は平均値によって値が変わってしまうため、が大きいほど分布が集中している、とは言えないという点である。
この点を改良した指数として、指数（Morisita's indexとも）がある。この指数はと同様にポアソン分布で1、一様分布で1より小、集中分布で1より大となるが、平均値に影響されず*2、平均値の異なる集団同士の集中度の比較に用いる事ができる。

ここでは番目の区画における個体数、は総個体数である。

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

参考文献

• 動物生態学研究法〈上巻〉 (1977年)

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

*2:一部例外はある

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

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

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

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

回帰式を変形させる

という回帰式を、

という風に変形する。このとき、は目的とするの値（上記例で言えば30）であるとすると、のときにとなるので、が求めるべきの値に等しいとういことになる。あとは、とについて、非線形回帰を用いて求め、の信頼区間を求める。
（参考）

• 逆推定: Pharmacometrics
• 最小2乗法，最尤法 線形モデル，非線形モデル

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

非線形回帰は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といういかにもな名前の関数が入っているが、それではイマイチ上手く行かなかった。使い方が間違っているのかもしれない。

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

事の発端

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

調べたこと

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

• 変数が空文字列かどうかを比較して確認する際は変数名をダブルクォーテーションで囲まないと止まる（全ては時の中に… : 【コマンドプロンプト】パラメータが空文字かどうか判定する）。
• コマンドの実行結果を変数に代入するには遅延環境変数（setlocal enabledelayedexpansionで使用可能になる!で囲まれた変数）for文を使ってゴニョゴニョする必要がある（バッチファイルで環境変数にコマンドの実行結果を格納(複数行) - 西海岸より）。
• echoで空行を作るにはecho;。
• ウィンドウを閉じないためにはpauseで一時処理を停止させる。

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

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

（※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)

グラフの横幅はずれない。
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)

ずれないが、下のグラフの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)

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)

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

追記

@nozma この件、rbind.gtable使うと楽ですよ。 grid.draw(rbind(g1, g2, size="first"))

— kohske (@kohske) 2015, 6月 2

ということを教えていただいた。
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の折りたたみ機能

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

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

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

上記4つは自動的に折りたたみ設定がされる。
例えば、関数や条件定義領域などのブレースで囲まれた部分があると、行数の横に三角が現れる。

これをクリックすることで折りたたみができる。

展開する場合は再び三角をクリックするか、折りたたまれている事を示している青いアイコンをクリックする。
（キーボードショートカットで操作する方法は用意されていないか、うまく動かないようだ。Ctrl + Command + L（Mac）で閉じられるときもあるし、閉じられないときもある。このコマンドは任意領域の折りたたみの際には上手く動作する。）

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

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

コードセクション

Rスクリプト中に、コメント記号（#）で開始し、4つ以上の-または=または#で終了する行があると、セクションの開始とみなされる。

どのセクション記号も機能は同様であり、セクション開始から次のセクション開始までが折りたたみの対象となる（＝サブセクションを作ることはできない）。
セクションは折りたたみの対象となるだけでなく、関数定義と同様にジャンプ機能の対象となり、その部分へ容易に移動できるようになる。

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

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を終了しても）が、ファイルを閉じて再度開いた際には全ての折りたたみが展開された状態となる。