BLOG

function関数の使い方

function関数

Rで複数の関数を組み合わせて値を算出しなければならないとき、function関数で定義してしまうと楽です。しかし、いったい何をやってるのかわかりにくいところもあります。ここでは、function関数を使った関数の定義について簡単に紹介します。

変動係数

変動係数(CV)とは、標準偏差を平均で除したもので、変動性の比較などでよく使われる指標です。RのデフォルトでCVは定義されていないのですが、比較的よく使うものなので、例として挙げてみます。まず、function関数の仕組みですが、引数の指定と、その引数を使った計算式の指定、という形になっています。

CV <- function(x){ sd(x)/mean(x)}

#example
> x <- rnorm(100,10,10)
> CV(x)
[1] 1.042615

上記のfunction()のかっこ内に引数(今回はx)として何を使うのか指定し、{}内に引数を使った計算式を定義します。

ベクター中の0要素をカウントする関数

少し込み入ったものとして、ベクター中に0がいくつかあるのかカウントし、そして0が何番目の要素かを返す関数を定義してみます。

zerocount <- function(x){
 n_zero <- length(which(x==0))
 id_zero <- which(x==0)
 return(list(n_zero=n_zero, id_zero=id_zero))
}

#example
>x <- rpois(100, 2)
> zerocount(x)
$n_zero
[1] 17

$id_zero
 [1] 12 16 21 22 25 31 35 40 43 67 74 77 78 79 86 87 95
 

こちらでは、返す要素が二つあるので、返り値をリスト化して返すよう指定しています(4行目)。

応用

発展的には、乱数発生などのプロセスを組み合わせることで、簡単なシミュレーションモデルを作ることができます。そのほか、自作の作図関数なども作ることができます。興味のある人はぜひトライしてみてください。

SPONSOR LINK

最頻値の推定関数

最頻値

平均、中央値、分散など、確率変数を特徴づける指標はいくつもあるけれども、多くのものについてはRで標準実装されている。しかし、地味に知りたい「最頻値」なるものは実装されていない。

最頻値の推定の仕方

離散値の場合は以下のコードで簡単に推定できる。

x <- rpois(100,10)
names(which.max(table(x)))

しかし、サンプルサイズが比較的小さい場合や、確率変数が連続値の場合にはとてもいい方法とは思えない。そこで、table関数の代わりに、density関数をつかった関数を定義する。

mfv <- function(x,from=min(x,na.rm=T),to=max(x,na.rm = T)){
  if(is.numeric(x)==F)stop("use numeric")
  tmp <- density(x,from=from,to=to)
  Mode <- tmp$x[which.max(tmp$y)]
  return(Mode)
}

# example
> x <- rnorm(1000,10,20)
> mfv(x)
[1] 12.16744

density関数は、確率変数の確率密度(のようなもの;実際はカーネル)を変数の傾度にそって推定してくれるもの。上で定義したものは、このカーネルの値が最も高くなる確率変数xの値を返す、という関数になる。細かいけど、いちいち計算するのは面倒なので、関数化してしまうのが楽。。。ちなみに、fromとtoに値を指定すると、その範囲の中だけでカーネルを推定する。例えば、カウントのような正の値しかとらない変数の最頻値を知りたいときは、from=0とすることで、正の実数の範囲の中だけで最頻値を計算してくれる。また、確率など0-1しかとらない、とわかっているのであれば、from=0、to=1とする。

だが、そもそもdensity関数がなにをどうやって推定しているのかようわかっていないので、正確ではない可能性がある。確認次第追記します。

SPONSOR LINK

「役立てるだけの能力がなかった」

ベイズの定理

統計の基本定理の一つである「ベイズの定理」が発表されたのは1763年。いまではGoogle検索のアルゴリズムにも採用されるほど重宝されているが(検索履歴を事前分布として使っているものと思われる;詳細は不明)、実際に実用化が進んだのは早くとも2000年以降だと思われる。ベイズの定理が発表された当時は、膨大な計算量が必要とされることから、「役に立たない」ものとみなされていたらしい。ベイズの定理の後に、Fisherにより最尤法なる画期的な推定手法が発表されると、その簡便性と正確性の高さから広い支持を得た。長い長い間、ベイズの定理が表立って利用されることはほとんどなかった。

マルコフ連鎖モンテカルロ

しかし、ごく最近になって、マルコフ連鎖モンテカルロなる逐次的な計算手法が開発された。これは、ベイズの定理と至極相性がよい。というのも、ベイズの定理における計算量の問題が、この計算手法によってかなり緩和されるからだ(※)。もちろん、手計算で行える範疇ではないが、近年のPCの計算能力をもってすれば、十分計算できるほどのものになる。「役に立たない」と思われていたベイズの定理に、やっと人類が追いついた。計算量の問題さえ切抜けてしまえば、ベイズのその他もろもろの利点が浮き彫りになり、各種研究分野で実用化が進んでいる。Google検索におけるベイズ推定の応用もそのひとつになるのだろう。

「今」役に立つ

「役に立つ」という言葉は、今目に見えている問題を相手にとった言葉なのだろう。役に立つ研究というものが、今後起こりうる未知の問題に対し、どれだけの示唆を与えてくれるのだろうか、と思う(「役に立つ研究」は、今の問題には役に立つけど、それっきりなんだろうなぁ)。「質の高い、けど今は役に立たない基礎研究」をサポートする姿勢というは、将来起こりうる不測の事態に対する先行投資とも思うのだけれど。税金を使っている限り、それは過ぎたお願いだろうか。参照:「なんの役に立つんですか?の暴力性

※ ベイズの定理の計算には、パラメーターpが任意の値をとるすべての事象の数え上げが必要になる。このため、普通の計算手法では、PCをもってしても不可能なほどの計算が必要になる。

ナダルが全仏で10回優勝したすごさについて

ナダルがテニスの全仏オープンで10回目の優勝を達成しました。すごすぎて、このすごさを表現するのに適切な言葉が見つかりません。なので、確率で表してみたいと思います。

指数分布のFitting

全仏オープンといえば、テニスの4大大会のひとつに数えられる大会。出場する選手は全力でここにスケジュールを調整してくるし、基本的にトップ選手はみんなエントリーする。なので、世界で一番優勝の難しい大会の一つといえる。ナダルはその大会で10回目の優勝を達成した。これが全人未踏なのはもちろんなのだが、どれくらいの確率で起きうるのかが気になる。なんとか推定できないか。そこでウィキペディアをググると、なんと全仏オープンの歴代優勝者があるではないか。優勝経験者のみを対象に、ナダルの偉業がどれほど「偉業」なのか推定してみる。横軸にある個人の優勝回数をとり、縦軸にそれぞれの優勝回数を達成した人数をとる(ヒストグラム)。これに対して、指数分布をあてはめる;

f(x) = θexp(-θx)
(xは優勝回数)

どれくらいの確率でおきるのか

指数分布でFittingすると、θ = 0.55と推定され、以下のような図になる。縦軸は確率密度で、値が高いほど対応するx軸の事象が生じる可能性が高い。

よく聞く名前では、フェデラー、ジョコビッチは1回、ヴィランデルは3回、ボルグが6回優勝。ナダルの「10回」という優勝回数が生じる確率は、0.004、つまり0.4%だった。上の図でいうと、(ほとんど見えないけど)赤で塗りつぶされているところの面積と対応する。ちなみにこれは、「全仏優勝経験者」が10回優勝する確率。なので、1000年全仏オープンを続けたら、確率論的には4人10回優勝する人がでるかもね、っていうくらい(実際にはこれよりもずっと低いと思う;選手生命の長さなども考えれば)。地殻変動レベルの現象。これはすごい。ナダル、おめでとう。

負の二項分布モデルの精度検証:glmer.nb、glmmamb、JAGS

精度はいかほどですか

この前紹介したけれども、lme4パッケージの更新で、ランダム効果入りの負の二項分布モデルが扱えるようになっていた。しかし、他の関数(glmmadmbやJAGS)と比べてどの程度の精度なのか、気になるところ。そこで、簡単に推定の精度比較をしてみる。モデルは極めて単純な線形モデル。

Y ~ NB(pi, n)
pi = n/(mui + n)
log(mui) = a + bXi + εj
εj ~ Normal(0,σR2)

まず、上記を真のモデルとするシミュレーションデータを生成し、それをglmer.nbなどに推定させる、という作業を100回繰り返す。100回分の推定値の頻度分布と真の値を比べるということをした。パラメーターの真の値はそれぞれa=1.0, b=0.5, σR = 0.7とした。

結果

赤い曲線が100回分の推定値の頻度分布。縦の赤い実線は100回分の推定値の中央値、グレーの実線は真の値。左から切片(a)、回帰係数(b)、ランダム効果のばらつき(σR)。

glmer.nb

glmmadmb

JAGS

こうしてみると、推定精度の観点から一番難がありそうなのはglmmadmb。特にランダム効果のばらつきの推定値は、右側の裾が随分と重い。これはランダム効果が二つ以上になったら絶望的ではなかろうか。glmer.nbは比較的安定していたが、回帰係数の推定値のピークが平らになっているのがやや不安。中央値は真の値に近いけれども。JAGSは可もなく不可もなく、推定としては一番マシかな、という気がする。ランダム効果が増えたときが気になるけれども(たぶんJAGSの推定がかなりまマシなものになる)、面倒なのでやらない。でも、最尤推定の負の二項分布モデルで、ランダム効果の入れ過ぎはよくないということは間違いなさそう。確認してないけど。

シミュレーションデータ生成に使ったコード

glmer.nbの精度検証に使ったコード。そのほかは割愛。
apply系を使えばfor構文が要らなくなる気がするけど。。。

for(i in 1:100){
  #simulate 100 data points from 10 groups with one continuous fixed effect variable
  x<-runif(100,0,1)
  f<-gl(n = 10,k = 10)
  data<-data.frame(x=x,f=f)
  modmat<-model.matrix(~x,data)
  
  #the fixed effect coefficient
  p <- NULL
  p[1] <- 1; p[2] <- 0.5
  fixed<-c(p[1],p[2])
  
  #the random effect
  p[3] <- 0.7
  rnd<-rnorm(10,0,p[3])
  
  #the simulated response values
  size <- 3
  mu <- exp(modmat%*%fixed+rnd[f])
  prob <- size/(mu+size)
  data$y<-rnbinom(100,size=size,prob=prob)

  fit <- glmer.nb(y~x+(1|f),data)
  est[[1]][i] <- fit@beta[1]
  est[[2]][i] <- fit@beta[2]
  est[[3]][i] <- fit@theta
}

スーパー堤防

父の実家が陸前高田なので、節目で被災地に訪れている(特にボランティアというわけではないです)。そのときに聞いた興味深い話。

被災以降、津波の「減災」を目的としたスーパー堤防がガンガン作られており、それはもうすごい勢い。生態学の観点からそれが好ましくないことはもちろんそうなのだが、現地に生きる人の感情を考えると極めて難しい問題だなぁ、と感じていた。しかし、それは大きな勘違い(みんな堤防ほしいだろう、という先入観)で、現地からでる「堤防」に対する意見は住民の間でも大きく割れており、そしてそのパターンがどうにも奇妙だった。

奇妙だったというのは、海の近くに住んでいる住人(つまり被災リスクが高い)ほど堤防反対の意見をもつことが多かったこと。逆に、海から離れた場所に住む人(被災リスクは相対的に低いと思われる)ほど堤防建設を熱望しているという。

海際の人は年配の人が多く住んでいる可能性あるので、年齢層との偽相関は拭えないが、海際の住人の共通した意見として、海に住んでいるのだからそういったリスクは承知の上、海を見れないことのほうがよっぽどつらいし、仕事もしにくい、とのこと(もちろん、これが一般的な傾向というつもりはございません。あくまで聞いた範疇です)。

まったく定量的な調査ではないですが、現地の人に聞く限り、上述のような傾向があるそうです。だれかきちんと証明しているのかな、これ。

さまざまな情報量基準

生態学の界隈では、情報量基準なるものが大人気です。そのくせ、(筆者も含め)それぞれの情報量基準なるものがどんな意味をもつのか、よくわかっていない場合がほとんどです。わかっている範疇で整理。

AIC(赤池情報量基準)

情報量基準といったらAICと言わんばかりですが、なんだったら決定係数との違いすらよく分かっていないかったのが博士課程のころの僕です。AICが低いモデルでは、「同じ母集団から同じサンプル数の別のデータセットを取り出してきたときにも、それなりのパフォーマンス(当てはまり)を発揮する」ということになる。だが、これはサンプル数を無限大に増やしたとしても、必ずしも「真のモデル」には収束しないらしい。あくまで評価基準が予測精度だからだろうか。

AICc

サンプルサイズが小さいときはAICc!みたいな使われ方をしている。だが実は、正規分布を仮定した補正をしているので、正規分布以外のモデルに対して適用するのはご法度らしい。Burnhamさんの本を参照のこと。

BIC(ベイズ情報量基準)

AICと異なり、こちらは真のモデルを選ぶらしい(あくまでサンプルが無限大になったとき、ということだろうか)。あまり使ったことはないが、たしかにAICとは全く異なるモデルを選ぶ傾向がある。予測精度ではなく、真に重要な変数を選び出したい場合にはBICを使った方がよいのだろう。式的には、AICの罰則項にサンプルサイズが入っていないのに対し、BICには入っている。そのため、サンプルサイズが大きくなったときに、前者は対数尤度の項の影響が相対的に大きくなり、より複雑なモデルを選びがちになるようだ。

DIC(Deviance Information Criteria)

こちらはベイズモデルなどでよく使われる情報量基準。だが、広く使われている割に、その数学的根拠は随分と希薄らしい(論文でこういった表記をよく見かける)。DICの式の構造としては、モデルの当てはまりに関わる指標(Deviance)に対し、パラメーター数で罰則を与えている。AICでは2k(kはパラメーター数)という罰則項があるが、これには数学的根拠がある。DICの場合は、そういった裏付けがない、ということだろうか。時間のあるときに調べよう。
また、DIC利用の際の注意点として、パラメーターの正規性が担保された状態でしか役に立たないらしい。詳細はよくわからないが、ベイズのご利益がでる階層線形モデルや状態空間モデルでは、この「パラメーターの正規性」がほとんど担保されないらしい。パラメーターの正規性とはなんぞや。。。(おそらく、パラメーターの事後分布が正規分布で近似できる、ということか)。だとしたらDICを使える場面などそうそうないのでは。

wAIC(widely applicable information criteria)

最近考案されたらしいwAICなるものは、階層をもった複雑なモデルにおいても適用できるらしい。これまたAICに引き続き日本人が発見したので驚き。日本人って統計につおいのだろうか。開発者のワタナベさんのページはこちら(ワタナベさんが開発したからwAICかと思った。違うのかorz)。内容をまったく理解していないのですが、浅はかな理解では予測精度を問題にしているような気がするので、モデル選択の結果の解釈としてはAICと同様になるような気がする。階層モデルのようなものにも適用できる点が違いか(あくまでユーザー視点からすれば)。もう少し勉強する。

lme4パッケージの更新

しばらく見ないうちに、最新版のlme4パッケージでは劇的な進化を遂げていたようなので、それについてメモ。

glmer.nbの出現

これまで負の二項分布を仮定したモデルで、ランダム効果をいれるとしたら、glmmADMBパッケージか、もしくはベイズモデル化(WinBUGSなど)する必要があった(ADMBもベイズみたいなもん?中身はよくわからない)。しかし、最新版のlme4パッケージでは、glmer.nbという関数が表れたことにより、負の二項分布でもランダム効果を簡単に入れることができるようになった。仕様については基本的なlmerと変わらない模様。ただし、頻度主義的に扱うには尤度関数が複雑すぎる気がするので、いったいどこまで結果が信用できるかは不明。>3のランダム項をいれるのはきっとご法度。

mcmcsamp()の廃止

昔のlme4には、最尤法ではなく、MCMCでベイズっぽく推定するための関数が用意されていた。しかし、最新版のlme4パッケージでは除外された模様。Rのhelpには以下のように書いてある。
Due to difficulty in constructing a version of mcmcsamp that was reliable even in cases where the estimated random effect variances were near zero (e.g. https://stat.ethz.ch/pipermail/r-sig-mixed-models/2009q4/003115.html), mcmcsamp has been withdrawn (or more precisely, not updated to work with lme4 versions >=1.0.0)

どうやらランダム効果のばらつきがゼロ付近のときに、あまりいい推定をしてくれないということで、やめにしたということらしい。個人的には、ランダム効果が多いとき(>3)の最尤推定は不安なので、mcmcsamp()の取り下げは手痛い、、、JAGSとかはコード書くのが面倒だからなぁ、、、

lmerにおけるREMLのデフォルト化

REMLとはなんぞや、と調べてみると、REstricted Maximum Likelihoodの略らしい。最新版のlmerでは、MLではなくREMLがデフォルトになっている。REMLとMLの違いは、(筆者の浅はかな理解によると)、前者は固定効果による自由度減少を考慮したうえでランダム効果の分散を推定し、後者は考慮しないらしい。つまり、後者の場合、ランダム効果のばらつきの推定は偏りが生じているということになる(不偏推定量でない?)。ただ、REMLではモデル選択ができないため(なぜかはよくわからない)、MuMInパッケージなどでモデル選択する場合はREML=Fとする必要がある。正規分布以外ではMLしか使えないので心配の必要なし。

glmerとlmerの(よくわからない)分離

正規分布の場合はlmer、それ以外はglmerということでよさそう。

QGIS: 属性テーブルの一括編集

QGISで属性テーブルの値を一括編集したいときの処理

1.変更対象の地物を選択する
2.属性テーブルからフィールド演算をひらく
3.既存カラムを更新を選択
4.変更したいカラムを選択
5.演算スペースに一括入力したい文字や数値を入力。
※文字の場合、シングルコーテーション(’ ‘)で囲む
※選択フィーチャーだけに適用のチェックボックスにチェック

ATMの悪夢

住み慣れた国を離れると、これまでの常識が通じないことも多々あり、いろいろと苦労しております。2017年の4月からアメリカに移動し、所属がUniversity of Minnesotaとなったのですが、さっそく新生活の立ち上げでひやり(というか実際損した)とした経験です。

日本ではお金を振り込んだり、預け入れをする際に、ATMの機械がその場で紙幣を数えてくれるのが普通だと思います。しかし、アメリカではそうではありませんでした。Depositのボタンを押すと、なにやら封筒がいるか、と聞いてくるのです。こちらの常識で考えると、当然封筒はいらないので、いらないと押し、そのまま紙幣を挿入口へ。ただウイーンといいながら吸い込まれていく紙幣。まったく数える気配がありません。よくわからないので、そのままその場を後にしたのですが、そもそも紙幣を数えていなかったこと、そして挿入金額が自己申告(途中で入力させられる)だったことに違和感を覚えるべきでした。不安になり、しばらく経ってからググってみると、、、げっ、、、封筒に金額と名前を記入して入れるだと、、、?オンラインアカウントの口座残高をみると、入金されていなことになっている(ちなみに500ドル)。なんとアメリカでは、ATMの預け入れは夜間金庫みたいなものになっており、行員があとでATMからお金を運び出し、手作業で数えているとのこと。なんじゃそりゃ。最近では日本のようなATMが大多数のようですが、運悪く私の使ったATMは旧式のもので、従来の方式が採用されていた模様。

そんな馬鹿な、と思い、支店にいって事情を話すと、行員の方が丁寧に対応してくれました。しかし時すでに遅し、紙幣をそのまま挿入したためか、500ドル入金したはずが480ドルになってしまいました。とほほ。全部なくなるよりはマシか。。。