Excelで複数条件抽出&集計

Excel

f:id:Rion778:20160713213905p:plain
このようなデータがあって、条件1と条件2に基づいて何らかの集計をしたいとする。

AVERAGEIFS

平均値を計算したいのであれば、AVERAGEIFS関数があるので、例えばこのように入力する。

=AVERAGEIFS(平均対象範囲,条件範囲1,条件1,条件範囲2,条件2)

最初の例で言えば、条件がこのように入力されているとして、
f:id:Rion778:20160713215530p:plain
入力は次のようになる。

=AVERAGEIFS(C:C,A:A,F3,B:B,G3)

AVERAGEIFSは高速に動作するので、列全体を範囲指定してもあまり問題が無い。
同様の集計関数にはCOUNTIFS、SUMIFSなどがあり、いずれも高速である。なお、Excel2007以降で使える。

配列数式

平均、カウント、合計以外の集計をしたい場合、悪名高き配列数式を使うという手段がある。
配列数式は配列を渡して配列が返ってくる仕組みで、一人で使っている分には便利だが、他人にファイルを引き継ごうとすると配列数式の説明を漏れなく行う必要があり、大変面倒くさいという代物である。
AVERAGEやVAR.Sなどのセル範囲を引数として受け取る関数は、セル範囲の代わりに配列を渡しても良いので、IFを組み合わせて配列を作成し、それを渡せば良い。

=AVERAGE(IF(IF(A1:A10000=F3,TRUE,FALSE)*
            IF(B1:B10000=G3,TRUE,FALSE),
            C1:C10000,""))

IF文中のTRUE、FALSEは1、0でも良い。
この式をCtrl+Shiftを押下しながら確定すると、配列数式に変換されて式の両端に{}が付く。見た目には{}以外の変化が無いので、大変分かりにくい。
そして配列数式の最大の欠点として、「遅い」というものがある。上記例ではわざわざセル範囲を1万行に制限してあるが、これをしないと1秒くらいの計算時間がかかってしまう。配列数式を使う場合、範囲指定は狭ければ狭いほど良い。

名前付き範囲

この配列数式の欠点を補う方法として、名前付き範囲を使用するという方法がある。
Excelにおける名前とは変数名のようなもので、指定のセル範囲を指定の名前で表記できるというものである。例えば、前述の式を次のように記述できる。

=AVERAGE(IF(IF(条件1=F3,1,0)*
            IF(条件2=G3,1,0),
            値,""))

これで大分可読性が向上する。
名前はリボンの[数式]から[名前の定義]で定義する。
ここで、次のように入力する。
f:id:Rion778:20160713222641p:plain
このようにOFFSET関数で参照範囲を定義することで、参照範囲が入力状況に対応して自動的に伸縮し、必要最低限の参照範囲となり、計算時間を節約できる。
なお、COUNTAで参照する列は、定義から想像出来ると思うが、全ての行に値が入っているような列でなければならない。MATCH関数等と組み合わせることで空白を許容することも可能であるが、長くなるのでここでは説明しない*1。
なお、名前には拡大倍率を40%以下にするとどこに名前を設定しているか一目で確認できるという便利機能があるが、上記の方法で設定した名前にはこの手段が使えない。
f:id:Rion778:20160713223803p:plain

まとめ

Rでやるべき。

*1:詳しい方法についてはExcel Hacks 第2版― プロが教える究極のテクニック140選に記述があるが、多分適当にググっても出てくるだろう。

グループ変数に応じて複数の折れ線グラフを書く

R

緑色の本のp.119に載っている図と似たものを書こうとして、

d <- data.frame(q = rep(c("q0.1", "q0.3", "q0.8"), c(9, 9, 9)),
                p = c(dbinom(0:8, 8, 0.1), dbinom(0:8, 8, 0.3), dbinom(0:8, 8, 0.8)),
                y = rep(0:8, 3))

というようなデータを準備したものの、baseパッケージではどうプロットしたら良いものか若干悩んだ。

baseの例

# base
matplot(reshape(d, timevar = "q", idvar = "y", direction = "wide")[, -1], 
        ylab = "p", xlab = "y", xaxt = "n", type = "b", pch = 16, lty = 1)
axis(1, at = 1:9, labels = 0:8)
legend("topright", inset = 0.05,
       legend = c("q = 0.1", "q = 0.3", "q = 0.8"),
       col = 1:3,
       pch = 1,
       box.lty = 0)

f:id:Rion778:20160705221323j:plain
linesで重ねる方法もあるけどいずれにせよ面倒だし使い回ししにくい。もっと良い方法はないものか。

続いてlatticeで描いてみた。

# lattice
library(lattice)
xyplot(p ~ y, data = d, group = q, type = "b", 
       auto.key = list(text =c("q = 0.1", "q = 0.3", "q = 0.8")))

f:id:Rion778:20160705221128j:plain
凡例の編集をしようと思わなければ書式が素直で一番理解しやすい様に思う。

ggplot2

# ggplot2
# install.packages("ggplot2")
library(ggplot2)
ggplot(data = d, aes(x = y, y = p, col = q)) + geom_point() + geom_line() +
  scale_color_hue(label = c(q0.1 = "q = 0.1", q0.3 = "q = 0.3", q0.8 = "q = 0.8"))

f:id:Rion778:20160705221214j:plain
latticeとさほど手間は変わらないけど、ggplot2はどうも覚えにくい気がする。初心者にaes()とかどう教えたら良いのかと毎回悩む。

RStudioのキーバインド変更

R RStudio AutoHotkey

Windowsでのお話。

Tools -> Modify Keyboard Shortcuts...

で出来るようになっているという話だけど、Ctrl+hをBackSpaceにしたいとかCtrl+zをPgUpにしたいとか思ってもどうもうまく出来ないので結局AutoHotkeyを使うことにした。ほんの少しゴニョゴニョしたいだけならAutoHotkeyが一番手軽な気がする。

  1. AutoHotkeyを入手&インストール(AutoHotkey)
  2. 拡張子.ahkなファイルに下記のスクリプトを記述して保存後ダブルクリック
#IfWinActive ahk_class Qt5QWindowIcon
^h::Send {BS}
^z::Send {PgUp}
#IfWinActive

ahk_classはAutoHotkeyに同梱のActive Window Info (Window Spy) で調べる。

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

R 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))

f:id:Rion778:20160111202035p:plain
凡例のラベルを操作する簡単な方法は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の指定が必要

f:id:Rion778:20160111212059p:plain

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

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

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

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

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

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

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

f:id:Rion778:20160111213151p:plain
labs()を使ってラベルだけを変更する場合も同様で、

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

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

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

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

facetのstripのtextを変更する

R

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

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

f:id:Rion778:20160103192059p:plain
このラベル部分は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)

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

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

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

飽和水蒸気圧の求め方

Meteorology R

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

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

飽和水蒸気量 - Wikipedia

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

Tetensの式はTetensの式ではない

Tetensの式やMurrayの式と呼ばれているこれらの式は、Magnusの式（後述するがこの呼称も適切ではない）と呼ばれている次式、
e_s=C \exp \left( \frac{A t}{B + t} \right)
のパラメータ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:

e_s(T)= 6.1094 \exp \left( \frac{17.625T}{T+243.04} \right) [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を含むいくつかの式の相対誤差をプロットすると、以下のようになる。
f:id:Rion778:20160103135736p:plain
日常的な気温の範囲ではどれも大差なく、むしろ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化＋その他いろいろ

Diary Windows

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

用意したもの

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インストール後の設定

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

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のブートローダーが立ち上がるという面倒な状態になってしまったのでこれを修正する。

bootrec /fixboot 
bootrec /fixmbr
exit

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

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

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