Ｒ備忘録 - 記事一覧
- 投稿者： みゅ
- カテゴリ： なし
- 優先度： 普通
- 状態： 完了
- 日時： 2009年06月30日 16時35分14秒
- 空間統計 ��
- 【gstat】バリオグラムから特定の位置における値を推定する
- 「The R Book」より引用
- 通常型クリギングシステムは
- という不偏性を確保した推定分散の最小化を解くことによって得られる以下のシステムを解くだけ、らしい．
- これで求まったｗをウエイトだと思って、観測されたデータを加重平均する．と、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])
Ｒ備忘録 /状態空間モデリング/donlp2/その他のメモ