sheephead 【統計】R

やっぱり駄目だった

１チェインで回り出したバージョンで、今度は３チェインで挑戦。
　がしかし、やはりtrap。もうどんな初期値でもだめ。iterationのグラフ見てみると、それぞれのチェインが素っ頓狂に飛んでいますな～。これはあかん。
事前分布もガッチガチなので、この時点で、もう完全に打つ手なし。

　こうなると、駄目駄目な因子突っ込んでも駄目って事で、一時潔く退くしかあるまい。あと、ぶっ込めそうな要因はあるが、またしても、めんどくせな計算が待っています。しかも、空間っぽい計算。

　しかし、ここまでイライラさせるArcGISなんて使いたくない。

　というわけで、とりあえずラスタをテキストで吐き出して、Ｒで計算させることにしました。これならイライラしても悪いのは、ほぼ自分なので精神衛生上すこぶる良い。
　

　というわけで、しばらくはwinbugsともお別れして、テキストとにらめっこの日々です。

2008年09月19日 | 【統計】R | トラックバック：0 | コメント：1

回り始めたzipクリギング

半ばあきらめぎみだったzipクリギングですが、タウの初期値諸々をいぢってやったらどうにか回り始めました。

　ただ、シカ多い少ないモデルに入ってるのは標高と傾斜と定数項、シカいるいないモデルに入ってるのは定数のみという超中途半端仕様で、しかもチェイン1こだけ。しかも事前分布ガッチガチ。
　とゆー問題ありありな感じですが、とりあえず、リモート妖精さんにチェイン数１でiteration長くしたバージョンの計算をお願いしてみる。

　んで、メイン妖精さんにはなかなか終わりそうにないArcGISでのモザイク処理。てか、今日気づいたんですがね、どうやらこのソフト、モザイクかけるファイルは同じフォルダにないと駄目っぽいですな。もーなんつーか、ビューワ的機能以外、このソフト引っかからずにまともに使えた試しがないんだが。

　もう、よっぽどRの方が賢い。何の文句も言わずDBFファイルも読み込めるし、shpも無問題。しかも、リモート妖精さんを使い始めて気づいたんですが、ネットワーク上のファイルもワーキングディレクトリにできるとゆー至れり尽くせりな仕様のワケで。

　ちなみに、lan内にsetwd()するときは以下のようにすればよいようです。

　setwd("//コンピュータ名/hogehoge/")
　ただ、linuxとかだと、ターミナル上でmountしないといけないとか、いろいろやり方も変わるそうな(ref.RjpWiki)

　２バイト文字も、たいがい大丈夫。いや、ほんとにすごいソフトだ。

　てわけで、現在使役できる妖精さんはフル稼働状態。ただ、腐れソフトにメイン妖精さん取られるのは、もったいなさ過ぎるよな。いちおデュアルマシンって事で。BUGSも2個動かせるし。

　そろそろ、GIS用マシンをこっそり作るか、サブ妖精さんを格上げするか・・・
　

2008年09月17日 | 【統計】R | トラックバック：0 | コメント：0

ベイジアンクリギングのつづきのつづき

ベイジアンクリギングのつづきのつづき
　
　とりあえず、昨日回していた計算は無事終わっていました。初期値も相当適当だったのに・・・
　学会の時は、trapの呪いにかかり、相当あくどいこともやったものですが、びっくりするくらいサクッと終了。ま、前やったのは、パラメータも結構多かったんですけど、今回はそんなにパラメータもないので、あまり迷うこともなかったのでしょう。

　収束判定やらなんやらは、とりあえず置いといて、glm()なアプローチ、つまり空間的自己相関を無視したやつと比べてみる。
　こっちは簡単でこんな感じ。

　　res <- glm(deer~elv,family=poisson)
　　res
　

　まずは、ベイジアンクリギングのsummary()

　　　　 2.5% 25% 50% 75% 97.5%
alpha　 1.174e-01 0.837025 1.160e+00 1.583e+00 2.457e+00
beta　 -2.256e-03 -0.001132 -4.664e-04 1.622e-04 9.949e-04
tau 　 1.389e-01 0.185900 2.175e-01 2.471e-01 3.144e-01
phi 　 5.943e-01 5.192250 1.049e+01 1.523e+01 1.956e+01
sigma 3.181e+00 4.047250 4.596e+00 5.380e+00 7.200e+00

　次はglmオブジェクト

Coefficients:
(Intercept) elv
2.9423393 -0.0009634

　ベイズジアンクリギングバージョンでパラメータが下がってるのは、空間的自己相関で説明してしまう分があるということなんでしょう。
　
　じゃあ、どっちがいいモデルなの？ってなると・・・
　気持ちとしては、何かしらの情報量基準でサクッと出したいとこではありますが、多分駄目な気がする。

　なら、交差確認法で、となるんですけどglm()あたりは、サクッと組めそうな気がしますが、MCMC計算入ってくるとそうはいかないわけで。
　なぜなら、R2winbugsでwinbugsはブートできるんでけど、終わらせることは多分できんのです。
　ということは、MCMC計算終わったとしたら、また手動でR2winbugsを実行させないといけなくなるという、ﾏﾝﾄﾞｸｾな状態になるはず。

　対策として考えられるのは、winbugsがマルチスレッドに対応していないというのと複数起動が可能なことを利用して並列処理させるという方法。（ODさんはそうしているそうな）
　ただ、これもクアッドコアとか、オクタコアとかレベルのCPUで旨味が出る話で、デュオクラスだとちとお話しにならない。

　ギークな方法としては、クラスターPCでやるという方法もあるみたいですけど（L. Tierney氏のsnowパッケージでクラスタ計算を行う - RjpWiki参照のこと）、これは個人的な力量不足と、職場のマシンをそういう使い方すると怒られそう、というので無理。

　てことで、可能性としてはチコチコMCMCさせては終了させ、実行、というプログラミングする人が最も嫌うルーチンな仕事をさせないといけなくなりそう。
　む～試行回数を何回にするかという問題はあれ、これは相当面倒くさい。

　なんか良い方法があるんだろうか。は～攻殻機動隊のデカトンケイルが欲しい・・・

2008年06月25日 | 【統計】R | トラックバック：0 | コメント：0

ベイジアンクリギング

新たなベイズ野望を果たすべく、新しく買った「Hierarchical Modeling and Analysis for Spatial Data」をちょー斜め読み。

　併せてkuboさんとこのぎょーむ日誌をマジ読みしてたら、ベイジアンクリギングなるワードがピキーンときてしまいました。どうやらWinbugsではSpatial.expなる関数で実現できそう。実はクラーク先生本にもベイジアンクリギングって書いてあるんだけど、込み入ったことはほとんど書いてなくて、もっぱらcar.normalにご執心なわけで。というわけで、「Hierarchical Modeling and Analysis for Spatial Data」、通称「矢印本」（表紙に矢印っぽいのがたくさん書いてある。たった今命名）を頑張って読んでみる。

　が、これまた非常に読みにくいんですな。大事なところになってくると数式がでてきたり、訳のわからない数学記号が出てきたりと相当ストレスが貯まります。なので、作者のサイトにあるBugsのコードを見ることに。うん、こっちの方がどう考えてもわかりやすい。

　http://www.biostat.umn.edu/~brad/data2.html

　これだと、spBayesとかより自由度が高く使えそうな気もします。

　参考になりそうなやつを探してみましたが、とりあえず日本語のやつは激少です。かろうじて、見つけたのが↓くらい。

　　空間統計モデルを用いたつくばエクスプレス沿線の地価の分析

　これが結構面白くて、距離のみに従属するベーシックなKrigingとbayesian Krigingを交差確認法で比較するというアプローチ。使ってるのもWinbugs。ただ交差確認法の試行回数は5回ってことで、ちょっと少ないんでね～の、な気もするもののMCMC計算にかかる時間も考えれば、妥当なとこかも。

　というわけで、ここ最近、RとWinbugs界隈をうろうろしながら、考えていた空間系ベイズの方向性もなんとなーく見えてきた感じ。

　てな話をODさんにメールしたら、car.normalな話題を10月頃発表するそうな・・・
　う～ん。あの悪魔的なベイズモデルの次にはもう空間系ですか。恐ろしい方です・・・

2008年06月21日 | 【統計】R | トラックバック：0 | コメント：0

九天社が倒産したらしい

Rの2chスレを見てたら、九天社が倒産したとの情報が。
　九天社と言えば、R関係で多数の書籍があって、僕も2冊ほど買わせていただきますた。
　む～、こうゆうニッチな出版社は、やっぱ波が激しいですな。

　欲しかったのに～って方は、在庫がなくなる前に急いで書店に注文されることをおすすめします。

　いちお、R関連の書籍の一覧を↓にかいときますけど、いずれなくなるかと。

　　http://www.9-ten.com/list/key_rtoukei.php

2008年06月14日 | 【統計】R | トラックバック：0 | コメント：0

統計読書の梅雨

梅雨入りして、調査もしばらくオフ。

　てことで、そんな時こそ貴重な自習の時間。てわけで、書籍まとめ買い。

　ベイズ関連2冊、R関連2冊、教養として1冊。

　必読はRプログラミングマニュアル。基本はRjpwikiのエキスを抽出したという感じですけど、改めて本で見ると結構知らないテクニックがあることに気づきます。スクリプトをバイトコンパイルして高速化とか。

　ベイズ関連では、5月末に出版された「マルコフ連鎖モンテカルロ法」はRのコードも盛りだくさんで、カエル本辺りの洋物を読まなくても、お手軽にベイジアンな手法がわかるというのは、ありがたいです。ただ、階層ベイズとか隠れマルコフとかcar.normalあたりのテクニックはほとんど触れてなくて、これまでの統計学をベイズ流にやるとこんな感じって感じか。
　それでもBrugsの具体的使い方も詳しく書いてて、R使いな方は一撃でマストの1冊。そういえば、初めて知ったんですが、Brugs使う時ってopenbugsｲﾝｽｺしなくていいんすか。てっきりｲﾝｽｺしないと使えないと思ってたからｲﾝｽｺしちまったい。

　とどめの1冊が「Hierarchical Mdeling and Analysis for Spatial Data」　タイトルそんままですが、ベイズ空間系の呪文がわんさか。適度にbugsのコードあたりも書いてますけど、基本数式。行列が出てきた日には、なんか偏頭痛がしてきます。
　ただ、避けて通れないので、頑張って読んでます。1.2マソもするしな。でも、car.normalの用途ってとても限定的なんぢゃないかって気がしてきた。もっと汎用性あるやつってないんだろか。あ、MCMC++使えってのは、なしね。

2008年06月10日 | 【統計】R | トラックバック：0 | コメント：0

geobugs

winbugsのメーリングリストでmaps2winbugsなるツールが紹介されてたので、早速ｲﾝｽｺ。

　SourceForge.net: maps2WinBUGS

う～ん。maps2winbugsというだけあって、いろんなベクターデータをwinbugsに読み込ませるために使うみたいですな。shpなんかも、お手軽にArcinfoなんかに変換してくれます。あと、バッファリングとかもできるみたい。つか、winbugsでベクターデータ表示させることができるなんて知らんかったですよ～
手持ちで、ほどよいサイズのベクターデータがなかったのでとりあえず、行政界のデータで練習してみました。
したら、びっくりするくらいサクッと、でけちゃったり。ちなみに林小班データ突っ込んだら固まっちゃいました。でかいデータはあんま良くないのかも・・・

http://blog-imgs-40.fc2.com/m/y/u/myuhe/bugs.png

動画付きのデモが付いてるので、僕のようなオバカさんでもサクッと使えます。

　んで、ついでにgeobugsのマニュアルをチラッと読んでみてたのですが、離散的なポイントデータでもcar.normalでいけそうな気がしてきました。

　要は、（多分）空間上のどことどこが隣接（近接？）しているかをはっきりさせれば良いわけですよね。

　どれくらいまでの距離を近接と言うんかい！！って怒られそうな気もしますけど・・・
　spbayesだとハイパーパラメータとかは設定できそうにないし・・・む～もうちょっとGeoBugsは勉強してみる価値あるかも。

2008年05月09日 | 【統計】R | トラックバック：0 | コメント：0

脱ESRI化計画

　次なる野望を果たすため、spBayes、maptoolあたりのお勉強。んでもって、関連書籍の注文をしてみたり。それにしても、ここんところのRの進化は異常というか、これまで人に説明するときは、タダの統計解析ソフトって説明してたんですが、地理情報関連のクラスもあるし、RGLなんか使った３次元表現もできて、WinbugsやGRASSも奴隷のように使役できるとなってくると、もう統計解析ソフトの範疇ではないような気が・・・
　データ解析環境統合プラットフォームとでもいいましょうか。

　とりあえず、ArcGISなる糞ソフトからの脱却を謀らねばなるまい。今日も保守関連の書類作っていたのですが、このソフト、信じられないバグだらけにも関わらず保守だけで年間20ﾏｿ以上も取りやがります。

　おいおい。保守する前にそのバグ、どげんかせえよ。詐欺に近いっつの。

2008年05月01日 | 【統計】R | トラックバック：0 | コメント：2

ベイズの苦悩

　最近、学会の時に作った階層ベイズモデルをこの期に及んで、いじろうかと努力しているところです。
　ま、初心者に毛の生えたような僕が突貫工事で仕上げたブツ故、問題ありありなのですが、そこを弱い頭をひねってどうにかしようとしているわけです。

　というわけで、とりあえず現時点での問題点を備忘録としてピックアップ。

・林床光環境問題
　地形による被陰と樹冠による被陰を別個に変数として入れても、そりゃあてはまり悪いのは当たり前。
　モデルの中で、どうにかして改良してやる必要あり。隠れ変数的にしてやるのもいいかもしれん。つっても、イマイチ隠れ変数のことわかっていないので、今後お勉強せんといかんです。

・時間スケール問題
　ま、林齢が林床植物群集を規定していると言うのは簡単なものの、林齢を固定効果な説明変数として使うのは良くないような気が。
　光環境を規定しているっていうのはあるが、それよりかは林床植物の人生（種子供給とか種間競争）みたいなものの影響が大きいのでは。
　ってすると、スケールを空間スケールから時空間スケールに拡張する必要があるかもしれん。つっても現時点では1時点データしかないので、2時点目取ったときのお楽しみか。でも1時点でもどうにかしたい気も・・・

・種間競争問題
　今回は、個体をざくっとポイントのような扱いにしてしまっているわけですが、当然個体は高さや種特性も違う。
　てことは、機能タイプや高さデータを能動的に入れてあげないといけないのかも。つっても、こればかりはノーアイデア。

　なことをぼ～んやり考えてたらkuboさんとこでも同じような話が出てたので、勢いで抜粋。

ふーむ …… 今までいいかげんにデータ解析されてた分野，とくに群集生態学なんかは今後は Bayes 化が進行していくのはマチガイないところだろうとおもうんだけど，このデータ解析でも見られるようになかなか面倒な部分がある．それは群集内の種間の相互作用の推定というやつだ．

生態学でいう相互作用とやらは (相関ではなく) 因果関係であり，因果関係を推定したいならば原理的には時間変化データが必要とされる．しかしながら，なかなかそのようにはデータがとられていない．データがないのに相互作用らしきモノを無理矢理に推定しようとするとこういう悪戦苦闘状態におちいる．

それじゃー気あいと根性で時間変化データとりゃあいいのか，というとハナシはそう簡単ではない．昆虫とかの動物群集の場合，群集の変化はわりと短期間に観察されるんだけど，季節変化とかイヤな外力に振りまわされる部分が大きいので，データとる回数が増えるにつれ次第に様相が複雑怪奇化していく，つまり自滅的な状況におちいったりすることなんかも．

いっぽうで植物群集なんかだと変化はもっとゆっくりしてるんだけど，ゆっくりしすぎていていったい何年ぶんのデータとれば言うところの相互作用とやらが推定されるの? ということになりがち．

ベイズ統計学の推定計算にもちいられる Markov chain Monte Carlo (MCMC) 計算，その故郷である統計力学でかくのごとき計算方法が発達してきた背景には「相互作用の強さの推定が困難 (だけど推定したい)」を克服することにあった．すなわち，システムの構成要素間のミクロな相互作用が推定困難な (何が原因で何が結果なのか観測による識別がしんどい) 状況のもとで，システムの内部状態の可能なパターンを列挙 &amp;amp; システム全体のマクロな統計量を算出できるところがウリになっている．

まあ，放棄スキー場植生の問題にたちかえると，ここで推定がうまくいってない理由は …… たとえばある quadrat で稚樹の本数がゼロだったときに，このモデルに組みこまれているふたつの説明
「ここは草本がたくさんいるから樹木稚樹は少ない」 (fixed effects 的)
「この場所にはそもそも樹木なんていないんだよ」 (random effects 的)
推定計算の中でこれらの「オレが正しい (オレが尤度を改善してやってるんだ)」バトルの決着がまったくさだかではなく，いつまでたっても「やっぱり草本による shading」 - 「いやいやそもそもいないんだ」の両者のあいだを激しくふらふらさまよっている (つまり熱力学的平衡とゆー範囲をこえて挙動に落ちつきがない) ことにありそうだ．

　てことで、自分の場合に立ち戻ってみると、GLMの延長上にランダム効果と固定効果のかけ算な要因を持ってくるのは、あんま良くないということで、良いんだろうか。ま、測って固定効果にしてしまえば良いんだろうが。あ、そういやシカが食べちゃって生えないということもあるんだった。ま、それはウンコを数えるので固定効果にしてしまえばいいか。

　なことを、ぐるぐる考えていくと、測ってないもんは結局わからないんだ、という極当たり前な結論に落ち着くわけです。だからといって、測りまくるというのも、ブルーカラーなお仕事の無限ループとなってしまう恐れが・・・

　はあ、妖精さんがそゆうめんどくさいことを全部してくれればいいのに・・・

2008年04月27日 | 【統計】R | トラックバック：0 | コメント：0

ポアソン分布の机上統計実験

　最近、統計の本とか読むようになり、改めて感じますのが、自分に数学的センスが皆無だということ。それとそのセンスって努力云々だけじぢゃ埋められないものもあるということ。
でも、見たりすればなんとな～く理解することはできるので、どうやって自分に視覚的に理解させるのかは、重要なポイントだったりします。

　というわけで、kuboさん的統計実験をやってみますた。お題は最近よく使うポアソン分布の過分散。それをばらつかせたら、どうばらつくのかな～という実験。スクリプトは↓な感じです。

test.x &amp;amp;amp;lt;- seq(0,990,by=10)
test.x &amp;amp;amp;lt;- rep(test.x,10)
id &amp;amp;amp;lt;- 1:1000
id &amp;amp;amp;lt;- as.character(id)
katamuki &amp;amp;amp;lt;- NULL
seppenn &amp;amp;amp;lt;- NULL

for(i in 1:100){

#test.y &amp;lt;- exp(0.005*test.x+1)
test.y &amp;lt;- exp(0.005*test.x+rnorm(length(test.x),1,0.3))
random.y &amp;lt;- rpois(length(test.x),lambda=test.y)
par(new=T)
plot(test.x,random.y,xlim=c(0,1000),ylim=c(0,500),xlab="標高",ylab="個体数密度",type="n")
res.lm &amp;lt;- glm(random.y~test.x,family=poisson)
res.coef &amp;lt;- coef(res.lm)
res.coef &amp;lt;- as.vector(res.coef)
katamuki &amp;lt;- append(katamuki,res.coef[2])
seppenn &amp;lt;- append(seppenn,res.coef[1])

par(new=T)
line &amp;lt;- function(x) {exp(res.coef[2]*x+res.coef[1])}
x &amp;lt;- seq(0,1000,length=100)
plot(x,line(x),col="red",xlim=c(0,1000),ylim=c(0,500),type="n",xlab="",ylab="")
lines(x,line(x),lwd=0.001,col="blue")

}