2011年11月30日

CPIIとCPIIIの間でのヒゲクジラの分布の変化

今日で11月終わるなぁと思って目標達成したかどうかを書こうとしたら、なんと、11月は目標立て忘れていた!確かに立て込んでいたけど・・・。来月はきちんと立てる。

Changes in baleen whale distribution patterns between CPII and CPIII. HIROTO MURASE, SHIGETOSHI NISHIWAKI, KOJI MATSUOKA, TAKASHI HAKAMADA AND TOSHIHIDE KITAKADO. (2011). SC/63/IA11.


国際捕鯨委員会(IWC)の科学小委員会(SC)に提出された論文。CPIIとCPIIIというのは、南極海で行っている南大洋鯨類生態総合調査計画の何週目に当たるのかを表す略語。それぞれ南極海の二週目と三週目を意味しているということ。それでこの論文はCPIIとCPIIIでクロミンククジラとその他のヒゲクジラ(ミナミセミクジラ、ナガスクジラ、ザトウクジラ)の分布がどのように変わったかを調べている。

使用した共変量は、CPIIかCPIIIかというダミー変数、水深、緯度と経度の二次の交互作用、そして水深800mからの距離と表面水温5度の等温線からの距離。調べているエリアは南極海を6分したときの番号である、Area I, IV, V-W。

まずArea Iでは、5度の等温線からの距離がCPIIとCPIIIで変わったことが分布に影響していることを示唆している。一般的に、夏の間クロミンククジラは氷縁に沿って分布することが良く知られている。それにも関らず、CPIIでは5度の等温線に沿って南極半島の沖合に多く分布していた。この論文では、5度の等温線を寒帯前線(Polar front)の指標として使っていた。そして南極半島はナンキョクオキアミの再生産するための場所としても知られている。つまり、前線から出来た渦が南極半島沖のナンキョクオキアミをより高密度に集めたために、そこにクロミンククジラも集まった可能性がある、としている。一方でCPIIIのときには前線が南極半島の沖合まで下がってこなかったために、同じ5度の等温線でも発見がなかった。そのため、CPIIとCPIIIでクロミンククジラの資源量に差があるのは、このような環境の違いも考慮して解析を行わなければならないであろう、としていた。そういう環境の違いや調査デザインがCPIIとCPIIIで違うことから、エリアごとの違いを考慮した解析のアプローチを考慮すべきだと結論づけている。クロミンククジラ以外の大型ヒゲクジラも資源量も増加し、分布も増加していることが確認されていた。

図が大変わかりやすく、クロミンククジラのいるグリッドセルの存在確率と、その他の鯨類がいる存在確率の差をとっている。そのおかげで差が大きい場所は住み分けているのではないか、と思えるような図になっている。また、CPIIよりもCPIIIのほうがクロミンククジラ以外の鯨類のいるセルが多くなっている。これは真似しようと思った。
ラベル: IWC 競争
posted by しばきん at 21:59| 神奈川 ☁| Comment(0) | TrackBack(0) | 研究生活 | このブログの読者になる | 更新情報をチェックする

2011年11月27日

誤差のある非線形生態動態システムのための統計的推論

Statistical inference for noisy nonlinear ecological dynamic systems. Simon N. wood. Nature vol466. 26. (2010).

僕はカオスのことはほとんどと言ってもいいくらい知らないので、勉強のために読んだ一本。といいつつ、ほとんどアイデアの勉強にしかならなかったのだけれど。
カオス的なふるまいをするある生物の個体数変動の時系列データが得られているとする。それをMCMCで最尤推定してやりましょう、という論文。カオスというと、なんだか混沌としたものを想像してしまうが、きちんとパラメタのある数式で記述出来る。ただし、超初期値依存そして予測不可能性を持つらしい。そんなもん、とうやって、最尤推定すんの、ということなのだが、そこを尤度の話に持って行ったこのWoodさんはすごい。要は自己回帰モデルを一回かますことで、自己回帰モデルのパラメタ間の比較の話にしてしまうのだ。つまり、

1.まず得られている生物の個体群動態モデルを決める。
2.そうすると、そこの動態モデルのパラメタが決まる(論文では、r, sigma, phi)
3.それらのパラメタを使う動態モデル「ではなく」、上記のデータに自己回帰モデルを当てはめる。
4.得られた自己回帰モデルのパラメタを真のパラメタと仮定する。
5.ここで、r, sigma, phiを勝手に決めて、個体群動態を表現する(観測データを勝手に作る)
7.勝手に作ったデータに同じ自己回帰モデルを当てはめる
8.これが、シミュレーションの一つ前の得られた自己回帰モデルのパラメタとどちらが先ほど仮定した真のパラメタと近いかを比べる。
9.7-8を繰り返してMCMC。
10.そうして得られたパラメタ群の平均値を生むような、r, sigma, phiが最尤推定値であろうよ、ということ。


このアイデアはすごい。
というわけで実際にやってみた。
簡単のため最小二乗法でやっていたり、負になるはずのないパラメタが負になるようなアルゴリズムになっているが、とりあえずそこそこ最尤推定出来るので無視(すんません)。以下、Core i5での実行時間。

##############
###真のパラメータ###
##############
set.seed(2)

###真のパラメタ###
r <- exp(3.8)
sigma <- 0.3
phi <- 10

###発生させる世代###
n <- 50
Nt <- numeric(50+1)

###観測誤差をポアソンで追加する前の真の動態###
Nt[1] <- 1
for (i in 1 : n) {
Nt[i+1] <- r*Nt[i]*exp(-Nt[i]+rnorm(1, 0, sigma))
}

###ポアソンで観測誤差###
Obs <- rpois((n+1), lambda=phi*Nt) #yt
plot(Obs, type="l")
points(Obs, pch=16)

###自己回帰モデルをかまして最尤推定出来るようにする###
###この自己回帰モデルのパラメータを推定することで###
###この後下で行うシミュレーションで得られる自己回帰パラメータと比較し###
###結果得られる最もここで推定する自己回帰パラメータに近かったr, sigma, phiが###
###最も尤もらしいという推論###
func_auto <- function (par, Obs)
{
beta1 <- par[1]
beta2 <- par[2]

###cubicスプラインを当てはめてパラメタ推定###
###最小二乗法###
return(
sum(
(
(Obs[2:(n+1)]^0.3) - (beta1*(Obs[1:n]^0.3) + beta2*(Obs[1:n]^0.6))
)^2
)
)
}

###上記のモデルをoptimで最尤推定###
res <- optim(func_auto, par=c(1, 10), Obs=Obs)
par_est <- res$par

#plot(Obs, type="l")
#points(Obs, pch=16)

#points(
# func_pred(par=par_est, Obs=Obs),
# col=2, type="l")

##############
###真のパラメータ###
##############

###MCMC###
k <- 500
par_est2 <- NULL
par_est2[[1]] <- c(1, 1)

###ステップ幅###
lag <- 0.25

###得られるパラメタのサンプルと初期値###
r_est <- sigma_est <- phi_est <- NULL
r_est[1] <- exp(3.8)-2
sigma_est[1] <- 3-1
phi_est[1] <- 10-2

###発生させるダミーの観測データ数###
Nr <- 500

system.time(
for (k in 1 : k) {

###パラメタの更新(本当は負の値にならないような工夫が必要)###
r_est[k+1] <- r_est[k] +runif(1, -lag, lag)
sigma_est[k+1] <- sigma_est[k] +runif(1, -lag, lag)
phi_est[k+1] <- phi_est[k] + runif(1, -lag, lag)

###以下上と同じ作業###
###r, sigma, phiのパラメタの候補を使って観測データを発生###
Obs2 <- NULL
par_2 <- matrix(NA, Nr, 2)
for(j in 1 : Nr) {
n <- 50
Nt2 <- numeric(50+1)

Nt2[1] <- 1
for (i in 1 : n) {
Nt2[i+1] <- r_est[k+1]*Nt2[i]*exp(-Nt2[i]+rnorm(1, 0, sigma_est[k+1]))
}

Obs2[[j]] <- rpois((n+1), lambda=phi_est[k+1]*Nt2) #yt

###同じ自己回帰モデルを当てはめる###
res <- optim(func_auto, par=c(1, 1), Obs=Obs2[[j]])
par_2[j, 1] <- res$par[1]
par_2[j, 2] <- res$par[2]
}

par_est2[[k+1]] <- c(median(par_2[, 1]), median(par_2[, 2]))

###メトロポリス法でパラメタの採択or棄却###
if(u[k] < sum((par_est-par_est2[[k+1]])^2)/
sum((par_est-par_est2[[k]])^2))
{
r_est[k+1] <- r_est[k+1]
sigma_est[k+1] <- sigma_est[k+1]
phi_est[k+1] <- phi_est[k+1]
} else {
r_est[k+1] <- r_est[k]
sigma_est[k+1] <- sigma_est[k]
phi_est[k+1] <- phi_est[k]
}
cat(k, "\n") #回数が不安なので今何回目のMCなのか確認
}
)
ユーザ システム 経過
607.61 1.18 612.43
###得られたパラメタの平均値###
> mean(r_est);mean(sigma_est);mean(phi_est)
[1] 44.56659
[1] 2.460081
[1] 9.372289

得られたパラメタのヒストグラム。
そこそこ尤もらしい山が見える。

MCMC_chaos_r.png

MCMC_chaos_sigma.png

MCMC_chaos_phi.png
posted by しばきん at 17:13| 神奈川 | Comment(0) | TrackBack(0) | 研究生活 | このブログの読者になる | 更新情報をチェックする

2011年11月20日

勇魚会雑感

昨日今日と松島で開かれた勇魚会に参加してきた。発表者には謝金が出ることを知らなかったので、思わぬ収入。もちろん、頂いた金額以上を宮城県で使ってきた。これも復興支援の一つだと考えている。

さて、残念ながら発表は賞を取ることが出来なかった。審査員のコメントを見ると、「数式がなかったのもあって大変分かりやすい発表であったが、それで結局その結果から言える結論はなんなのか」という主旨のコメントが多かった。こういうコメントをいただくのは狙い通りであったのだが、少し読み間違った。

そもそも勇魚会とは学生主体の会だと聞いていたので、聴衆もおそらく学生が多く、発表をするとしても、それはゼミの延長線上に当たるのだろうと踏んでいた(加えるなら、イルカ大好き!可愛い!キャー!みたいな女子学生が多く、個体数推定?なにそれ美味しいの?と言われると踏んでいた)。そのため、ゼミや個人的な指導のときに、学生に考えさせるような感じでスライドを作ってしまっていた。つまり、「あえて結論は載せない、しかし結論自体は仄めかしておいてあとは自分で考えさせる」、という具合である(ちなみにこの技は人から盗んだものです。上手く噛み合えば面白い議論が続いて大変楽しいです)。だから、学生が多いならそこから何が言えるのか一緒に考えましょうかー、といった具合に持っていく予定であった。しかし蓋を開けてみると、学生も多いが研究者の方々も多く、コメントにもあるようにこれは采配をミスったかな、と反省。そりゃまぁ、結論をあえて載せていないんだから、当然そこを疑問に思うわけです。本当なら、はいそこの箇所は確かに疑問ですよね、じゃあその結果って何を意味していると思いますか?、くらいから始めて一緒に議論することで、本当の意味で僕の言いたいことが伝わるはずだったが、完全に計算ミス。今思えば傲慢かつ阿呆な采配だった。

これが学会に準ずるものであるなら、やはり結論をちゃんと述べる必要があると思う。相手に仄めかすだけで考えさせる、というのはその場合傲慢以外の何物でもないであろうから。一緒に考えようぜー等とふるまうのは、ただ時間の無駄である。そういうわけでその点は次に生かすべき点であると思う。反省。まぁ僕が何をやっているか何が出来るのか、はわかりやすく伝わったはずなので、個人的には満足である。学生にも伝わったみたいだし。






まぁ結論としては、ランダムサンプリングとシミュレーションって大事ですね、ということです。これはこれではしょりすぎ。
posted by しばきん at 22:09| 神奈川 ☁| Comment(0) | TrackBack(0) | 研究生活 | このブログの読者になる | 更新情報をチェックする

2011年11月19日

分子系統樹推定の理論と実践

統計数理研究所のセミナー「分子系統樹推定の理論と実践」に参加してきた。自分の分野とは特に関係は今のところないと思い込んでいるが、勉強のため。ようやく、自分で系統樹の最尤推定が出来る気がしてきた。とりあえずの暫定プログラムを載せておく。ただし、まだ正の方向に10%のシステマティックなエラーがあるので、どこかにコードミスがあることが確認されている。講師の二人に聞いたが、バグ取りならず。どこが悪いのだろう。近日中に修正する。


unit <- c("ATCG")
n <- 200
human <- rep(unit, n)
human <- unlist(strsplit(human, split=""))
res2 <- NULL
N <- 1000
for (i in 1 : N){
chimp <- sample(c("A", "T", "C", "G"), length(human), rep=T)
human_chimp <- table(human, chimp)
sum(human==chimp)
sum(diag(human_chimp))
sum(human_chimp)-sum(diag(human_chimp))
ga <- table(c(human, chimp))[1]/sum(table(c(human, chimp)))
gc <- table(c(human, chimp))[2]/sum(table(c(human, chimp)))
gg <- table(c(human, chimp))[3]/sum(table(c(human, chimp)))
gt <- table(c(human, chimp))[4]/sum(table(c(human, chimp)))
#ga;gc;gg;gt
func4 <- function (x, y, v1) {
-1*
log(
(
sum(human==chimp)/sum(human_chimp)*
1/4 + (1-1/4)*exp(-v1)
)^x +
(
(sum(human_chimp)-sum(human==chimp))/sum(human_chimp)*
1/4*(1-exp(-v1))
)^y
)
}
res <- optimize(f=func4, c(1, 2, 3),
x = sum(diag(human_chimp)),
y = sum(human_chimp)-sum(diag(human_chimp)))
res2[i] <- res$objective
}
hist(res2)
length(human)
ラベル:統計
posted by しばきん at 09:57| 神奈川 ☔| Comment(0) | TrackBack(0) | 研究生活 | このブログの読者になる | 更新情報をチェックする

2011年11月15日

別冊図書館戦争I・II

別冊図書館戦争I


別冊図書館戦争II


2つまとめて、何年ぶりに読んだかであろうにやにやが止まらない本であった。論文書きの休憩時間に読むものではないな、と思った、良い意味で。この本のシリーズに柴崎という女性のキャラクターがいるのだが、どうも僕はこのキャラがいたく気に入ったらしい。もともと、読書した後、そのキャラに大きく影響を受けるタチなのだが(例えば村上春樹の作品を読んだ後、しばらく思考が作中のそれに近くなる)、今は柴崎に影響を受けている。普段苦手な、こう言ったら相手がこう考えるだろうな、等の相手を思い遣る思考に今日はえらい頭が働いた。おまけに、僕は一日に一回くらい将棋ゲームを休憩中にするのだが、おかげさまで(?)今日は全勝だった。とは言っても、そんな感じに影響されている状況は半日から一日で収まる。面白いのがそのときにする仕事は普段と少し違った出来になるような気がしていて、プラスに働くこともマイナスに働くこともある。新書や新聞を読んでもこうはならないので、要は無批判に本を読むとこうなるのかなぁ、とか思っている。
posted by しばきん at 23:45| 神奈川 ☁| Comment(0) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする
×

この広告は1年以上新しい記事の投稿がないブログに表示されております。