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

こないだに引き続き、ベイジアンクリギング自習の続き。
　日本語の資料はとりあえずググってみたので、google scholarで英語のやつを捜してみる・・・んですが、science directに契約していないので、とりあえずタダで見られるもの限定です。

　geobugsのマニュアル見てみると、hierarchicallyに書くやり方と、Non-hierarchicallyに書くやり方があってhierarchicallyに書くやり方をおすすめしているよう。
　つまり、クリギングの条件付き分散のパラメータをハイパーパラメータとして制御するってことのよう。
　う～ん、となると、相当自由度高く書けるという理解でいいのか。とすると、ますます、spatial.expって有用な気がしてきた。あえて強引に使いどころを考えてみると

　　spatialなデータ　-> spatial.exp
spatio-tempralなデータ　-> car.normal

なんですかねえ。点過程解析なんかと対比しても、やっぱcar.normalって特殊な気がするしなあ。

　と、ここまで進んで、とても重要なことに気づく。

　そういえば、クリギングって何っすか？

　・・・

　というわけで、3歩進んで、5歩ほど戻って、「空間データモデリング」を読んで、クリギングの勉強してみる。
　が、内容が難しく、理解できず。行列とか出てきた時点で、思考停止するし。

　早々と、あきらめて、次は先日倒産した九天社出版のRbookを見てみる。こちらは、解説が易しいのでなんとな～くわかってきました。

　ふむふむ、未知のポイントを既知のポイントから距離重み付けの値で推定する・・・んだと思う。

　ただですね、spatial.exp()で推定するパラメータとクリギングで推定するパラメータとの対応がイマイチわからんのですよ。

　ナゲットとかシルとかレンジとかは、結局わからんってことで良いんだろうか・・・

　ま、わからんことだらけなんですが、当座winbugsの書き方くらいはイメージできてきたので、早速書いてみます。

まずは、Rにはこんな感じで書く。

library(R2WinBUGS)
library(coda)

setwd(“d:/【研究】シカ関連/ベイジアンクリギング”)
dat<-read.csv(“2006生息密度調査.csv”,header=T)

deer <- as.numeric(dat$deer)
x <- as.numeric(dat$x)
y <- as.numeric(dat$y)
elv <- as.numeric(dat$elv)
n <- length(deer)

setwd(“d:/Rscript”)

data <- list(“deer”,”x”,”y”,”elv”,”n”)
parameters <- c(“alpha”,”beta”,”tau”,”phi”,”sigma”)
inits <- list()
inits[[1]] <- list(alpha=-13.78, beta=50, tau=0.01, phi=5.0 ,w=rep(1,times=n))
inits[[2]] <- list(alpha=-13.78, beta=100, tau=2.28, phi=4.9 ,w=rep(2,times=n))
inits[[3]] <- list(alpha=-13.78, beta=200, tau=0.001, phi=5.1 ,w=rep(3,times=n))
test <- bugs(data, inits, parameters, model.file=”krig.bug”, n.chains=3, n.iter=10000, n.burnin=8000, debug=T)

そして、Winbugsには、こんな感じに・・・

model {
# Likelihood
# Loop for Plot
for (i in 1:n) {
# Loop for Species
deer[i] ~ dpois(out[i])
log(out[i]) <- alpha + beta*elv[i] + w[i]
mu[i] <- 0
}

w[1:n] ~ spatial.exp(mu[],x[],y[],tau,phi,1)

# Parameters and hyper parameters

alpha ~ dnorm(0,0.001)
beta ~ dnorm(0,0.001)
tau ~ dgamma(0.001, 0.001)
phi ~ dunif(0.05, 20)

sigma <- 1/tau

}

　う～～ん、なんか相当適当ですが、bugsが、ガリガリ計算しているので多分大丈夫だろう。
　しばらく、見てるけれども計算は一向に終わりそうにないです・・・
　時間がかかりそうなので、結果は明日見てみようかしらん。