トップ 差分 一覧 ソース 検索 ヘルプ RSS ログイン

BugTrack-R備忘録/41

R備忘録 /状態空間モデリング/donlp2/その他のメモ

R備忘録 - 記事一覧

空間統計 - 【gstat】バリオグラムから特定の位置における値を推定する

  • 投稿者: みゅ
  • カテゴリ: なし
  • 優先度: 普通
  • 状態: 完了
  • 日時: 2009年06月30日 16時35分14秒

内容

  • 空間統計
  • 【gstat】バリオグラムから特定の位置における値を推定する

推定

  • 「The R Book」より引用
  • 通常型クリギングシステムは
  • という不偏性を確保した推定分散の最小化を解くことによって得られる以下のシステムを解くだけ、らしい.
  • これで求まったwをウエイトだと思って、観測されたデータを加重平均する.と、x0の位置の推定値が求まる.
library(gstat)
library(lattice)
# - メウス川
data(meuse)
data(meuse.grid)
x <- variogram(object=log(zinc)~1, locations=~x+y, data=meuse, cutoff=3000)
plot(x)
model1 <- vgm(psill=0.6, model="Sph", range=600, nugget=0)
plot(x, model=model1)
model2 <- fit.variogram(object=x, model=vgm(psill=0.6, model="Sph", range=600, nugget=0))
plot(x, model=model2)
g <- gstat(id="logzn", formula=log(zinc)~1, locations=~x+y, data=meuse, model=model2)
xo <- predict.gstat(g, meuse.grid[,1:2])
# - 推定
krig_mat <- diag(0, nrow(meuse))
krig_mat[lower.tri(krig_mat)] <- variogramLine(model2, dist_vector=dist(meuse[,c("x","y")]))[,2]
krig_mat <- t(krig_mat)
krig_mat[lower.tri(krig_mat)] <- variogramLine(model2, dist_vector=dist(meuse[,c("x","y")]))[,2]
krig_mat <- rbind(krig_mat, 1)
krig_mat <- cbind(krig_mat, 1)
krig_mat[nrow(meuse)+1,nrow(meuse)+1] <- 0

solve_krig_mat <- solve(krig_mat)
y <-
apply(meuse.grid[,c("x","y")], 1,
    function(x, meuse, solve_krig_mat){
        x <- c(x[1], x[2])
        n <- nrow(meuse)
        ret <- variogramLine(model2, dist_vector=sqrt((meuse[,"x"] - x[1])^2 + (meuse[,"y"] - x[2])^2))
        krig_vec <- c(ret[,2], 1)
        ret <- solve_krig_mat %*% krig_vec
        sum(log(meuse$zinc) * ret[1:n])
    },
    meuse=meuse, solve_krig_mat=solve_krig_mat
)
mean(xo[,3] - y)
plot(xo[,3], y)
  • このように一緒になった
  • まあgstatなら以下ようにすればいいだけなんだが.
g <- gstat(id="logzn", formula=log(zinc)~1, locations=~x+y, data=meuse, model=model2)
xo <- predict.gstat(g, meuse.grid[,1:2])

コメント