単変量の分布母数の推定からあれこれ:fitdistrplusパッケージが便利
ある確率分布を手持ちの単変量データに当てはめ,その母数を最尤推定するという方法があって,Rでは通常MASSパッケージのfitdistr関数を使うのがお手軽なんだけど,fitdistrplusという便利なパッケージが出ていてこれがとてもいいかも。ま,結局はoptimにかぶせているだけだが。
例えば,
- さまざまな可視化
- ブートストラップ
- 任意の母数を固定した状態での別の母数の推定
- 最尤推定以外の方法
- 適合度比較
みたいなのをかなり手軽にやってくれる。いちいち自分で書かなくてもいい。
可視化
たとえば,可視化に関してならば,
#数値例の作成 library(fitdistrplus) dat<-rgamma(300,2,3) #分布の確認 hist(dat) plot(ecdf(dat)) #これをガンマ分布に当てはめてみる fit<-fitdist(dat,"gamma") plot(fit) summary(fit)
ほら,便利。これだけ。
左上は,ヒストグラムに推定された母数の値を項にもつ確率密度関数。
右上はQ-Qプロット。この例だと,何故か大きい値があんま当てはまっていないことがわかる。
左下は経験累積密度関数と推定された母数の値の累積密度関数。右下はP-Pプロット。
もちろん,summaryも見やすい。こんな感じで普通欲しいのは一気に出してくれる。
Fitting of the distribution ' gamma ' by maximum likelihood Parameters : estimate Std. Error shape 1.965111 0.1488415 rate 2.841928 0.2450182 Loglikelihood: -156.0496 AIC: 316.0992 BIC: 323.5068 Correlation matrix: shape rate shape 1.0000000 0.8785199 rate 0.8785199 1.0000000
対数尤度プロットなる機能もある。これは,2つの母数の組み合わせの尤度をヒートマップみたいにしたのね。
尤度関数の2次元版,みたいなもんだろうか。ガンマは二母数だからね。
llplot(fit)
ブートストラップ
ブートストラップで誤差や信頼区間の検討もお手軽にできる。
fit<-fitdist(dat,"gamma") #B=1000でパラメトリック bfit1<-bootdist(fit,bootmethod="param",niter=1000) #B=500でノンパラメトリック bfit2<-bootdist(fit,bootmethod="nonparam",niter=500) #結果の表示 summary(bfit1) summary(bfit2)
ここではブートストラップした値の中央値と2.5%と97.5%点を出してくれる。←これはパーセンタイル法でのブートストラップ信頼区間でもある。
結果を色々見るには,
#パーセンタイル法での信頼区間 bfit1$CI #ブートストラップ推定値 boot.shape<-bfit1$estim[,1] boot.shape boot.rate<-bfit1$estim[,2] boot.rate #これをヒストグラムで描いて,中央値,パーセンタイル信頼区間を描き入れる。 par(mfrow=c(1,2)) hist(boot.shape,col="lightgray") abline(v=quantile(boot.shape,c(0.025,.5,.975)),col=2) hist(boot.rate,col="lightgray") abline(v=quantile(boot.rate,c(0.025,.5,.975)),col=2)
片方の母数を固定してもう片方の母数を最尤推定
やりたかったのこれ。ある理論的な見地から,片方の母数,例えば形状母数を任意の値に固定した状態で尺度母数を推定したいとか,これって結構外国語教育研究的によくあること。こういうときは,
#形状母数を無理やり3に固定(ホントは2) fitfix<-fitdist(dat,"gamma",method="mle",fix.arg=list("shape"=3)) #結果を見る summary(fitfix) Fitting of the distribution ' gamma ' by maximum likelihood Parameters : estimate Std. Error rate 4.525801 0.15086 Fixed parameters: value shape 3 Loglikelihood: -149.1229 AIC: 300.2459 BIC: 303.9496 #1つだから尺度母数のみの対数尤度曲線が出る llplot(fitfix)
最尤推定以外の推定方法
最尤推定以外の方法も実装されている。methodに入れてやればいい。
- mle 最尤推定。もちろん初期値も与えられるし,optimの方法も指定できる
- mme モーメント法
- qme 分位点法
- mge 最大適合度法 よく知らない
fitfist(dat,"gamma",method="mme") fitfist(dat,"gamma",method="qme") fitfist(dat,"gamma",method="mge")
適合度比較
件のデータを,ワイブル,対数正規分布にも当てはめてみて,適合度を比較する。
#フィットさせる gam<-fitdist(dat,"gamma") weib<-fitdist(dat,"gamma") logn<-fitdist(dat,"lnorm") #比較(リストで入れる) gofstat(list(gam, weib, logn), fitnames=c("gamma", "Weibull", "lognormal")) #結果 goodness-of-fit statistics gamma Weibull lognormal Kolmogorov-Smirnov statistic 0.04626154 0.04626154 0.08208341 Cramer-von Mises statistic 0.05906861 0.05906861 0.43467475 Anderson-Darling statistic 0.40584196 0.40584196 2.93641451 Goodness-of-fit criteria gamma Weibull lognormal Aikake's Information Criterion 278.3222 278.3222 321.8015 Bayesian Information Criterion 285.7298 285.7298 329.2090
こんな感じ。便利ね。コルモゴロフ・スミルノフやアンダーソン・ダーリング検定ってのはわかる。残りのクラメール・フォン=ミーゼスってのは,こんなこと?なるほど。
Cramér–von Mises criterion - Wikipedia
…ま,最初はこうやってモデリングについて慣れていくことが大事で,こういうのがまさに出発点って感じだなあと最近しみじみ。
最尤推定した母数のもとでの確率密度曲線をヒストグラムに描き足す
えっと,リクエストがあったのでここに書く。(SPSSならそれっぽい曲線もつけてくれるのにRはそんなこともできないのか?といわれた)
まずはこんなデータがあるとしよう。形状母数が3,尺度母数1のガンマ分布にしたがう300個。
set.seed(0) dat<-rgamma(300,3,1)
これをヒストグラムにする。
hist(dat,col="lightblue",main="")
まずはガンマ分布を仮定して,母数を最尤推定する。
library(MASS) fitdistr(dat,densfun="gamma")
ま,3.01と1.03とかが推定されたと。
この母数のもとでの関数をヒストグラムに描き足したい。
まずはhist関数に,freq=Fを入れる。上に少し余裕を持たせるためにylimを指定しよう。
hist(dat,col="lightblue",main="",freq=F,ylim=c(0,.4))
すると縦軸が密度になる。
この上に関数を描けばいい。
たとえば,
x<-seq(0,9,.001) lines(x,dgamma(x,3.01,1.03),col="blue")
はい。上手く当てはまっている感ある。
これで,別の関数も試してみたらいいかもだ。ワイブルと正規分布いこう。
#ここまでの作図 hist(dat,col="lightblue",main="",freq=F,ylim=c(0,.4)) x<-seq(0,9,.001) lines(x,dgamma(x,3.01,1.03),col="blue") #ワイブルからパラメーターを取る w.shape<-fitdistr(dat,densfun="weibull")$estimate[[1]] w.scale<-fitdistr(dat,densfun="weibull")$estimate[[2]] #ワイブルを描き足す lines(x,dweibull(x,w.shape,w.scale),col="red") #正規分布からパラメーターを取る mu<-fitdistr(dat,densfun="normal")$estimate[[1]] sigma<-fitdistr(dat,densfun="normal")$estimate[[2]] #正規分布を描き足す lines(x,dnorm(x,mu,sigma),col="green")
緑がダメ,みたいな感じがわかる。
ちなみにAICを比較してみたり…
aics<-numeric(3) aics[1]<-AIC(fitdistr(dat,densfun="gamma")) aics[2]<-AIC(fitdistr(dat,densfun="weibull")) aics[3]<-AIC(fitdistr(dat,densfun="normal")) aics
ほい。
野良観測変数?
なんていうんだろ,それが何を測るかは一切分かっていないけど,その変数の利活用を暗に強いられているような変数。
オンライン学習履歴データとかがまさにそうだ。オンライン教材のログイン時間とか。これは何か測定を目指す構成概念がこのデータに先んじて導出されているわけではない。従事率?学習者の能力?真面目さ?それを考える前に記録するよう実装されているし,なんていうか「何かに使うために実装しておけ,単純に考えて有益そうだが,それが何のためかは俺シラネ」感に満ちている。
そういう測定する構成概念とか,むしろその測定モデルに関する研究者の合意(ある意味妥当性の検証)がないまま,この観測変数を使うっきゃない,みたいな感じになっている。
およそデータマイニングといわれるような手法やそれを利活用する分野は大体そうだ。ある店に客が頻繁に来る確率を顧客ロイヤリティといったりするのは知っているけど,これは,顧客ロイヤリティという潜在変数を測定するために人工的に開発された変数ではない。店に来たり,高い買い物をしたりすることが先にある。
別にこれが悪いとか,そういう変数はダメだ,という話では全くなくて。データサイエンスってのはそういうものだし。
でも,こういう何を測るかわからない変数と,構成概念が先んじて理論的な見地より導出されて,測定モデルがしっかりあって人工的に作成された変数(優れたテストや質問紙や,広く受け入れられた実験手法の結果変数)は違うね,ってだけ。
いうならば,こういうのは野良観測変数。
それで,外国語教育研究も他の分野に漏れずに,その野良観測変数が最近一気に増えた。
…私はただ,野良猫を人々がさも自分のペットであるかのように勝手に好きな名前をつけて呼び合う,というありふれたことにならないかということが心配だ。
そして,ある人が「タマ」と呼ぶ猫を,別の人がそれを見て「え,これはミーちゃん」ってなって,そして「どちらかといえばタマじゃね?」「ミーちゃんだし」とか,そんなわけわからん話になって結局その猫はみんなのペットでもなく結局は野良猫のまんまで結局猫可愛がりということばのように愛情なるものも結構無責任だなってなってなんていうか猫は救われてなきゃだめなんだ。自由で,主体性を体現しているかのように。そして,バラの名前は,なんでも同じ香りがするって話もあるけどな。
オンラインでできるPOS Tagger
なんか最近こういうのやたら多いね。POSってはpart of speech。品詞のこと。
ウェブブラウザでテキスト投げてやればPOSをタグして返してくれるっていうお手軽サービス。
1. Parts-of-speech.Info - POS tagging online
2. NLTK POS Tagging - API & Demo | Text Analysis Online | TextAnalysis
4. https://cogcomp.cs.illinois.edu/page/demo_view/POS
5. Free English Part-of-Speech Tagging Service
6. LX-Center
分散共分散行列を相関係数行列に変換するR関数
自分で関数書けばいいってのもごもっともだし,他にも色々あるんだけど,bayesmパッケージにnmatっていう関数がある。
vcovm<-matrix(c(3,2,2,3),2,2) cm<-nmat(vcovm) matrix(cm,2,2)
手法上の利点とノウハウの違い
手法上の利点,または手法的に優れている(methodologically sound)手法というのは,基本的に適用する対象と独立してその有用性が与えられるものだ。だから,「手法上の」(methodologically)という言葉をつける。
これは主に数理的な,または技術的な意味で従来のものを包含してしまうか(俗にいう一般化された),または明らかにその手法自体のアフォーダンスに当たる要素について,客観的で,もはや自明ともいえるような評価が与えられるものに限る。
私は1980年代の生まれなので,ビデオデッキもしっているし,その規格争いについても知っている。そして新しいDVD規格,そしてBlueray規格もしっている。おそらく,これら新しいもの(記録方法など)は古いものに比べて手法的に優れている。保存できる情報量とか,そういう数値でそれが評価できる。
でも,大事なことは,我々がしばしば道具や手法のよしあしを議論するときに,「適用する対象と独立しているか」について見逃してしまうことだ。
「昭和のアダルトビデオを見るのはやはり擦り切れてフィルム焼けしたビデオに限るな,ブルーレイや動画配信なんて好きじゃない。これは俺の青春だ」っていうちょっと気持ちの悪いおっさんがどこか広島あたりにいてもいいし,スマフォよりもガラケーの方が好きなひともいていいし,誰が竹枕で寝てもいいし。…ああ,そうですか,となる。たいていの人は。
…これはノウハウだ。手法的に(技術的に)ブルーレイがビデオよりも優れているとか,そういうのとはちょっと関係の薄いことだ。
私の分野,外国語教育研究は計量的な分析について非常に後進的だが,「手法Aより手法Bの方がいい」みたいな話でいっぱいだ。結構なことだ。
でも,「手法的に優れている」ことを主張するためには,あくまでも「その手法を応用する対象」とは独立して,客観的な方法で示したらいい。たとえば統計モデルでは,従来のモデルを新しいモデルが含んでいることを数理的に説明できたらいい。
でも一方のノウハウは,ある対象における文脈の中の話だ。
分散分析より重回帰がいいというと話をよく聞く。大概のひとは耳だこだ。最近も聞いた。手法上の利点の話をすれば,正直どっちもあまり変わらない。成立過程を問わなければ,どちらも一般線形モデルの仲間と捉えられるし,一般線形モデルは一般化線形モデルへ一般化された。混合モデルをこれに適用できるようになったし,もはやMCMCなどの数値解析法に頼れば,従属変数,独立変数,階層レベルでの母数の分布,そういったものはとても自由にモデリングできるようになった。こういった発達を考えれば,手法上そんなに差はない。それに分散分析や重回帰は構造方程式モデリングで代用できるし。
手法上のことをいえば,典型的な分散分析は,独立変数がカテゴリカルで,典型的な重回帰よりも独立変数の数が少ない。…手法上の観点でいえば,でも,まあそれくらいだ。
この話の問題はノウハウに当たることだろう。
独立変数の尺度水準は,単純に研究対象による。
実験を計画し,条件を割り当てたりしたとき,その尺度水準がカテゴリカルなのは自明であるし,場合によってはカテゴリカルな尺度水準の方が理論上望ましい,または理論上それらしい時もある。観測変数とその背後にある変数のそれぞれの水準は,基本的には理論または必要性に関する議論から得られるべきだ。
また,独立変数の数についても,ランダムサンプリングができており,割り当てた水準との交絡が無い限り,従属変数に対する交絡変数の影響は無視できる。そういった時に無駄に変数を増やすことは研究計画の一貫性を著しく下げる。
もっと簡単にいうと,統計モデルを選択するときに,ある研究においてそれが分散分析モデルに帰着したからといって,それがなんだというんだろう。
こういったある目的下における親和性などを無視して,さも手法上の利点であるかのようにある手法をアピールすることは違う。そしてある具体的な状況下・目的を踏まえないノウハウなんてもはやなに。
統計の手法はより制約が自由に,そして理論的見地からモデルがより導出しやすく,そしてよりユーザーフレンドリーになる方法で発達していく。しかしその中で,例えばその古い方法の制約が満たされて,十分に理論的見地よりそれで検証できることが明白ならば,別に新しい手法使わなくていいんじゃない。そういうノウハウも重要でしょう。手法が自由になる,というのはそういうことでしょう。
新しい手法,新しい手法というようにとらわれるのも,もはや1つの理解不足だな,というように自分にも常々言い聞かせたいものだ。
もちろん,ある場合によって重回帰の方が断然いいだろうというノウハウには大いに賛同だ。そしてこれはもっとも議論する余地のないところだと思う。