2010年12月15日

自作通常クリギングプログラム

行列でパラメータ推定することを覚えたので、クリギングを自作プログラミングしてみました。なるほど、いつもクリギングは時間がかかるんだなぁと思ってましたが、こういう計算をしているのであれば納得です。

参考文献
クリギングとその地理的応用

> ####################
> ###自作クリギングプログラム###
> ####################
> #pdf("kriging.pdf")
> set.seed(2)
> library(lattice)
> library(spatstat)
> library(gstat)
> ###生データの作成###
> raw <- data.frame(x=rnorm(100, 0, 1), y=rnorm(100, 0, 1), z=0)
> raw[, 3] <- log(10*abs(1/(raw[, 1]*raw[, 2])))
> raw2 <- raw
>
> ###生データのサンプル場所の確認###
> #plot(raw2[, 1], raw2[, 2])
>
>
> ###まずは既存のパッケージを使ってやってみる###
>
> ###座標系としてやる###
> coordinates(raw2) = ~ x + y
>
> ###適当にあてはめ###
> v <- variogram(z~1, raw2)
> #plot(v)
>
> ###初期値を与えて当てはめ###
> v.fit <- fit.variogram(v, vgm(psill=0.9, "Gau", range=7, nugget=0.2))
>
> ###シルとかレンジの値###
> #summary(v.fit)
> v.fit
model psill range
1 Nug 0.1704907 0.0000000
2 Gau 1.8520082 0.6434817
>
> ###点に含まれる数が等しくないと当てはまりが悪いように見えることもある###
> #plot(v, model=v.fit)
>
> ###予測グリッドの作成###
> mesh.grid <- expand.grid(x=seq(-3, 3, 0.1), y=seq(-3, 3, 0.1))
>
> ###クリギング補間###
> gs <- gstat(id="z", formula=z~1, data=raw2, model=v.fit)
> coordinates(mesh.grid) = ~x+y
> res_gs <- predict.gstat(gs,mesh.grid)
[using ordinary kriging]
> aa <- data.frame(z=res_gs$z.pred ,expand.grid(x=seq(-3, 3, 0.1),
+ y=seq(-3, 3, 0.1)))
> ###推定された図###
> levelplot(aa[, 1] ~ aa[, 2] + aa[, 3], main="package_used")
>
> ###平均値と分散を確認###
> mean(aa[, 1])
[1] 2.661233
> var(aa[, 1])
[1] 0.7555209
>
>
> ###次に自作でプログラムを組んで結果が一緒になるか確認###
>
> ###地点間の距離と分散を計算###
> d_temp <- u_temp <- v_temp <- a <- a2 <- NULL
> for(i in 1 : nrow(raw)) {
+ d_temp <- sqrt((raw[i, 1]-raw[, 1])^2+(raw[i, 2]-raw[, 2])^2)
+ u_temp <- (raw[i, 3]+raw[, 3])/2
+ v_temp <- (raw[i, 3]-u_temp)^2+(raw[, 3]-u_temp)^2
+ a <- cbind(d_temp, v_temp)
+ a2 <- rbind(a2, a)
+ }
>
> ###バリオグラムにモデルを当てはめるときのカットポイントを指定###
> ###ここが違うので既存パッケージの結果とぴったり一緒のバリオグラムモデルは推定されないことに注意###
> bp <- seq(0, 2.1, length=16)
>
> ###切った間の平均値を計算してそれにモデルを当てはめる###
> a2_temp <- a3 <- a4 <- NULL
> for(i in 1 : (length(bp)-1)){
+ a2_temp <- subset(a2, a2[, 1]<=bp[i+1]&a2[, 1]>=bp[i])
+ a3 <- apply(a2_temp, 2, mean)
+ a4 <- rbind(a4, a3)
+ }
>
> ###確認###
> #plot(a4, pch=16)
>
> ###ガウスモデルのパラメータ推定###
> func_vario <- function(par, data, dis) {
+ inter <- par[1]
+ phi1 <- par[2]
+ phi2 <- par[3]
+ return(
+ sum(
+ (data-(inter + phi1*(1-exp(-(dis/phi2)^2))))^2
+ )
+ )
+ }
> res <- optim(c(1, 2, 1), func_vario, data=a4[, 2], dis=a4[, 1], hessian=T)
>
> ###出力用に関数を作る###
> func_vario_res <- function(inter, phi1, phi2, dis) {
+ return(
+ inter + phi1*(1-exp(-(dis/phi2)^2))
+ )
+ }
>
> ###さっきのプロットに線を重ねる###
> #points(a4[, 1], func_vario_res(inter=res$par[1], phi1=res$par[2],
> # phi2=res$par[3], dis=a4[, 1]), type="l", lwd=2
> #)
>
> ###任意の地点間の距離を計算してくれる関数の作成###
> func_dist <- function(pred, obs) {
+ pred_p <- ppp(pred[, 1], pred[, 2], c(-1e+5, 1e+5), c(-1e+5, 1e+5))
+ obs_p <- ppp(obs[, 1], obs[, 2], c(-1e+5, 1e+5), c(-1e+5, 1e+5))
+ return(crossdist(obs_p, pred_p))
+ }
>
> ###生データ間の距離を計算###
> dis <- func_dist(pred=raw, obs=raw)
>
> ###距離を使ってバリオグラムの計算###
> X <- func_vario_res(inter=res$par[1], phi1=res$par[2],
+ phi2=res$par[3], dis=dis)
>
> ###重みの合計は1###
> X <- rbind(X, 1)
> X <- cbind(X, 1)
>
> ###対角成分は0距離###
> diag(X) <- 0
>
> ###予測グリッドの作成###
> mesh.grid <- expand.grid(x=seq(-3, 3, 0.1), y=seq(-3, 3, 0.1))
>
> ###予測グリッドと生データの距離を計算###
> dis2 <- func_dist(pred=mesh.grid, obs=raw)
>
> ###バリオグラムの計算###
> for(i in 1 : nrow(mesh.grid)){
+ X2 <- func_vario_res(inter=res$par[1], phi1=res$par[2],
+ phi2=res$par[3], dis=dis2[, i])
+ X2 <- c(X2, 1)
+
+ ###クリギング方程式を解く###
+ ###一番下はラグランジュの未定係数法のパラメータであって推定には使わない###
+ mesh.grid[i, 3] <-sum(raw[, 3]*as.numeric(X2%*%solve(X))[1:100])
+ }
> #x11()
> levelplot(mesh.grid[, 3] ~ mesh.grid[, 1] + mesh.grid[, 2], main="original_program")
>
> ###パッケージで推定したものとほとんど同じ平均値と分散###
> mean(mesh.grid[, 3])
[1] 2.662301
> var(mesh.grid[, 3])
[1] 0.7765148

> ###4分位点も確認。大体一緒###
> quantile(aa[, 1])
0% 25% 50% 75% 100%
1.028423 2.156357 2.553124 2.760064 6.581620
> quantile(mesh.grid[, 3])
0% 25% 50% 75% 100%
0.974662 2.154865 2.559820 2.748471 6.761476
> #dev.off()
> ####################
> ###自作クリギングプログラム###
> ####################


出来た図。kriging.pdf
パッケージを使ったものとほとんど同じ結果。なぜ違うかというと、バリオグラムのパラメータ推定時のカットポイントの値が違うから(だと思う)。

posted by しばきん at 20:23| 神奈川 ☁| Comment(0) | TrackBack(0) | 空間統計 | このブログの読者になる | 更新情報をチェックする
×

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