クリギング入門ー空間データ推定の確率論的アプローチ　阪田義隆先生　コロナ社

```r
#p.74 理論バリオグラム(最小二乗法でのあてはめ)
#実際に当てはめのデータを作成する際には近しい区間で区切って
#値を平均化したものを使用する

data(EuStockMarkets)

u=EuStockMarkets[1:200,1:4]
hs=array(0,dim=c(nrow(u),nrow(u)))
for(i in 1:nrow(hs)){
hs[i,]=apply((t(u)-c(u[i,]))^2,2,sum)  
}

r=2+(3-2)*(1-exp(-hs/100000))
r=r*100

Model=function(x){
h=x[1:(nrow(r)^2)];x=x[-c(1:(nrow(r)^2))]
a=x[1];b=x[2];c=x[3]  
if(model==1){#Sherical Model
predicts=ifelse(h<a,b+(c-b)*(3*h/(2*a)-0.5*(h/a)^3),c)
}
if(model==2){#Exponential Model
predicts=b+(c-b)*(1-exp(-h/a))  
}
if(model==3){#Gaussian Model
predicts=b+(c-b)*(1-exp(-(h/a)^2))  
}
if(model==4){#Hole-Effect Cyclic Model
predicts=b+(c-b)*(1-cos(pi*h/a))
}
return(predicts)
}

f=function(x){
vec=c(c(hs),x)
predicts=Model(vec)  
return(predicts)  
}

model=1
Y=c(r)
h=0.001
lam=0.01
ite=10^4
sita=c(mean(hs[hs>0]),rep(0.01,2))
eta=10^(-2)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(Y),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*beta)
print(sum((Y-f(sita))^2))
}

```


```r

#4．クリギング
#単純クリギング　p.86

#各座標
x=c(2,1,4,6,5.5)
y=c(4,3,2.5,2,4)
u=cbind(x,y)

#観測ベクトル
z=c(5.3,4.5,4.6,2.9,3.2)

#距離
h=dist(u,method="euclidean")
hd=array(0,dim=c(length(z),length(z)))
hd[lower.tri(hd)]=c(h);hd=hd+t(hd)

#推定点１
h1=c(1.5,3)
h0=sqrt(apply(t(t(u)-h1)^2,1,sum))

#非類似度と距離　p.79
r=dist(z,method="euclidean")
rd=array(0,dim=c(length(z),length(z)))
rd[lower.tri(rd)]=c(r);rd=rd+t(rd)/2
plot(c(h),c(r),xlab="h",ylab="r")

m=mean(z)

#理論バリオグラム
#球関数モデル
rhk=1.84-1.84*(3*hd/(2*7.69)-(1/2)*(hd/7.69)^3)
rhk[hd>7.69]=0
rhk0=1.84-1.84*(3*h0/(2*7.69)-(1/2)*(h0/7.69)^3)
rhk0[h0>7.69]=0
wk=solve(rhk)%*%rhk0
ZSKk=m+sum((z-m)*wk)
sigSKk=1.84-sum(wk*rhk0)
pthik=c(-1,wk)
Ck=rbind(c(1.84,rhk0),cbind(rhk0,rhk))
costk=t(pthik)%*%Ck%*%t(t(pthik))


#ガウスモデル
rhg=1.64*(1-exp(-(Inf/2.91)^2))-1.64*(1-exp(-(hd/2.91)^2))
rhg0=1.64*(1-exp(-(Inf/2.91)^2))-1.64*(1-exp(-(h0/2.91)^2))
wg=solve(rhg)%*%rhg0
ZSKg=m+sum((z-m)*wg)
sigSKg=1.64*(1-exp(-(Inf/2.91)^2))-sum(wg*rhg0)
pthig=c(-1,wg)
Cg=rbind(c(1.64*(1-exp(-(Inf/2.91)^2)),rhg0),cbind(rhg0,rhg))
costg=t(pthig)%*%Cg%*%t(t(pthig))

#指数モデル
rhe=2.5*(1-exp(-Inf/5.53))-2.5*(1-exp(-hd/5.53))
rhe0=2.5*(1-exp(-Inf/5.53))-2.5*(1-exp(-h0/5.53))
we=solve(rhe)%*%rhe0
ZSKe=m+sum((z-m)*we)
sigSKe=2.5*(1-exp(-Inf/5.53))-sum(we*rhe0)
pthie=c(-1,we)
Ce=rbind(c(2.5*(1-exp(-Inf/5.53)),rhe0),cbind(rhe0,rhe))
coste=t(pthie)%*%Ce%*%t(t(pthie))



```

```r

#通常クリギング

#各座標
x=c(2,1,4,6,5.5)
y=c(4,3,2.5,2,4)
u=cbind(x,y)

#観測ベクトル
z=c(5.3,4.5,4.6,2.9,3.2)

#距離
h=dist(u,method="euclidean")
hd=array(0,dim=c(length(z),length(z)))
hd[lower.tri(hd)]=c(h);hd=hd+t(hd)

#推定点１
h1=c(1.5,3)
h0=sqrt(apply(t(t(u)-h1)^2,1,sum))

#非類似度と距離　p.79
r=dist(z,method="euclidean")
rd=array(0,dim=c(length(z),length(z)))
rd[lower.tri(rd)]=c(r);rd=rd+t(rd)/2
plot(c(h),c(r),xlab="h",ylab="r")

m=mean(z)

#理論バリオグラム
#球関数モデル
rhk=1.84*(3*hd/(2*7.69)-(1/2)*(hd/7.69)^3)
rhk[hd>7.69]=1.84
rhk=cbind(rbind(rhk,rep(1,ncol(rhk))),c(rep(1,nrow(rhk)),0))
rhk0=1.84*(3*h0/(2*7.69)-(1/2)*(h0/7.69)^3)
rhk0[h0>7.69]=1.84
wk=solve(rhk)%*%c(rhk0,1)
muk=wk[length(wk)]
wk=wk[-length(wk)]
ZSKk=sum(z*wk)
sigSKk=muk+sum(wk*rhk0)


#ガウスモデル
rhg=1.64*(1-exp(-(hd/2.91)^2))
rhg=cbind(rbind(rhg,rep(1,ncol(rhg))),c(rep(1,nrow(rhg)),0))
rhg0=1.64*(1-exp(-(h0/2.91)^2))
wg=solve(rhg)%*%c(rhg0,1)
mug=wg[length(wg)]
wg=wg[-length(wg)]
ZSKg=sum(z*wg)
sigSKg=mug+sum(wg*rhg0)


#指数モデル
rhe=2.5*(1-exp(-hd/5.53))
rhe=cbind(rbind(rhe,rep(1,ncol(rhe))),c(rep(1,nrow(rhe)),0))
rhe0=2.5*(1-exp(-h0/5.53))
we=solve(rhe)%*%c(rhe0,1)
mue=we[length(we)]
we=we[-length(we)]
ZSKe=sum(z*we)
sigSKe=mue+sum(we*rhe0)


```

```r


#平均値クリギング

#各座標
x=c(2,1,4,6,5.5)
y=c(4,3,2.5,2,4)
u=cbind(x,y)

#観測ベクトル
z=c(5.3,4.5,4.6,2.9,3.2)

#距離
h=dist(u,method="euclidean")
hd=array(0,dim=c(length(z),length(z)))
hd[lower.tri(hd)]=c(h);hd=hd+t(hd)

#推定点１
h1=c(1.5,3)
h0=sqrt(apply(t(t(u)-h1)^2,1,sum))

#非類似度と距離　p.79
r=dist(z,method="euclidean")
rd=array(0,dim=c(length(z),length(z)))
rd[lower.tri(rd)]=c(r);rd=rd+t(rd)/2
plot(c(h),c(r),xlab="h",ylab="r")

m=mean(z)

#理論バリオグラム
#球関数モデル
rhk=1.84*(3*hd/(2*7.69)-(1/2)*(hd/7.69)^3)
rhk[hd>7.69]=1.84
rhk=cbind(rbind(rhk,rep(1,ncol(rhk))),c(rep(-1,nrow(rhk)),0))
wk=solve(rhk)%*%c(rep(0,nrow(hd)),1)
muk=wk[length(wk)]
wk=wk[-length(wk)]
ZSKk=sum(z*wk)
sigSKk=muk


#ガウスモデル
rhg=1.64*(1-exp(-(hd/2.91)^2))
rhg=cbind(rbind(rhg,rep(1,ncol(rhg))),c(rep(-1,nrow(rhg)),0))
wg=solve(rhg)%*%c(rep(0,nrow(hd)),1)
mug=wg[length(wg)]
wg=wg[-length(wg)]
ZSKg=sum(z*wg)
sigSKg=mug


#指数モデル
rhe=2.5*(1-exp(-hd/5.53))
rhe=cbind(rbind(rhe,rep(1,ncol(rhe))),c(rep(-1,nrow(rhe)),0))
rhe0=2.5*(1-exp(-h0/5.53))
we=solve(rhe)%*%c(rep(0,nrow(hd)),1)
mue=we[length(we)]
we=we[-length(we)]
ZSKe=sum(z*we)
sigSKe=mue




```


```r

#非定常クリギング
#普偏クリギング

#各座標
x=c(2,1,4,6,5.5)
y=c(4,3,2.5,2,4)
u=cbind(x,y)
X=cbind(rep(1,length(x)),u)

#観測ベクトル
z=c(5.3,4.5,4.6,2.9,3.2)

#回帰パラメータ
beta=solve(t(X)%*%X)%*%t(X)%*%z
f=t(X)*c(beta)

#距離
h=dist(u,method="euclidean")
hd=array(0,dim=c(length(z),length(z)))
hd[lower.tri(hd)]=c(h);hd=hd+t(hd)

#推定点１
h1=c(1.5,3)
h0=sqrt(apply(t(t(u)-h1)^2,1,sum))
f0=c(1,h1)*beta

#非類似度と距離　p.79
r=dist(z,method="euclidean")
rd=array(0,dim=c(length(z),length(z)))
rd[lower.tri(rd)]=c(r);rd=rd+t(rd)/2
plot(c(h),c(r),xlab="h",ylab="r")

m=mean(z)

#理論バリオグラム
#球関数モデル
rhk=1.84-1.84*(3*hd/(2*7.69)-(1/2)*(hd/7.69)^3)
rhk[hd>7.69]=0
rhk=cbind(rbind(rhk,f),rbind(t(f),array(0,dim=c(ncol(X),ncol(X)))))
rhk0=1.84-1.84*(3*h0/(2*7.69)-(1/2)*(h0/7.69)^3)
rhk[h0>7.69]=0
wk=solve(rhk)%*%c(rhk0,f0)
muk=wk[-c(1:nrow(hd))]
wk=wk[c(1:nrow(hd))]
ZSKk=sum(z*wk)
sigSKk=1.84-sum(c(rhk0,f0)*c(wk,muk))


#ガウスモデル
rhg=1.64-1.64*(1-exp(-(hd/2.91)^2))
rhg=cbind(rbind(rhg,f),rbind(t(f),array(0,dim=c(ncol(X),ncol(X)))))
rhg0=1.64-1.64*(1-exp(-(h0/2.91)^2))
wg=solve(rhg)%*%c(rhg0,f0)
mug=wg[-c(1:nrow(hd))]
wg=wg[c(1:nrow(hd))]
ZSKg=sum(z*wg)
sigSKg=1.64-sum(c(rhg0,f0)*c(wg,mug))


#指数モデル
rhe=2.5-2.5*(1-exp(-hd/5.53))
rhe=cbind(rbind(rhe,f),rbind(t(f),array(0,dim=c(ncol(X),ncol(X)))))
rhe0=2.5-2.5*(1-exp(-h0/5.53))
we=solve(rhe)%*%c(rhe0,f0)
mue=we[-c(1:nrow(hd))]
we=we[c(1:nrow(hd))]
ZSKe=sum(z*we)
sigSKe=2.5-sum(c(rhe0,f0)*c(we,mue))

```

```r

#非定常クリギング
#外生ドリフトクリギング

#各座標
x=c(2,1,4,6,5.5)
y=c(4,3,2.5,2,4)
u=cbind(x,y)


#観測ベクトル
z=c(5.3,4.5,4.6,2.9,3.2)

#回帰パラメータ
s=c(3.85,4.17,3.22,2.65,2.77)
b=cov(s,z)/var(s)
a=mean(z)-mean(s)*b
S=t(cbind(rep(1,length(s)),s))
R2=cor(z,a+b*s)^2

#距離
h=dist(u,method="euclidean")
hd=array(0,dim=c(length(z),length(z)))
hd[lower.tri(hd)]=c(h);hd=hd+t(hd)

#推定点１
h1=c(1.5,3)
s1=c(1,4.02)
h0=sqrt(apply(t(t(u)-h1)^2,1,sum))




#理論バリオグラム
#球関数モデル
rhk=1.84-1.84*(3*hd/(2*7.69)-(1/2)*(hd/7.69)^3)
rhk[hd>7.69]=0
rhk=cbind(rbind(rhk,S),rbind(t(S),array(0,dim=c(nrow(S),nrow(S)))))
rhk0=1.84-1.84*(3*h0/(2*7.69)-(1/2)*(h0/7.69)^3)
rhk[h0>7.69]=0
wk=solve(rhk)%*%c(rhk0,s1)
muk=wk[-c(1:nrow(hd))]
wk=wk[c(1:nrow(hd))]
ZSKk=sum(z*wk)
sigSKk=1.84-sum(c(rhk0,s1)*c(wk,muk))


#ガウスモデル
rhg=1.64-1.64*(1-exp(-(hd/2.91)^2))
rhg=cbind(rbind(rhg,S),rbind(t(S),array(0,dim=c(nrow(S),nrow(S)))))
rhg0=1.64-1.64*(1-exp(-(h0/2.91)^2))
wg=solve(rhg)%*%c(rhg0,s1)
mug=wg[-c(1:nrow(hd))]
wg=wg[c(1:nrow(hd))]
ZSKg=sum(z*wg)
sigSKg=1.64-sum(c(rhg0,s1)*c(wg,mug))


#指数モデル
rhe=2.5-2.5*(1-exp(-hd/5.53))
rhe=cbind(rbind(rhe,S),rbind(t(S),array(0,dim=c(nrow(S),nrow(S)))))
rhe0=2.5-2.5*(1-exp(-h0/5.53))
we=solve(rhe)%*%c(rhe0,s1)
mue=we[-c(1:nrow(hd))]
we=we[c(1:nrow(hd))]
ZSKe=sum(z*we)
sigSKe=2.5-sum(c(rhe0,s1)*c(we,mue))

```

```r

#非定常クリギング
#コクリギング

#各座標
x=c(2,1,4,6,5.5)
y=c(4,3,2.5,2,4)
u1=cbind(x,y)
x=c(0.5,1,2,4)
y=c(4.5,1,2,4)
u2=cbind(x,y)

#観測ベクトル
z1=c(5.3,4.5,4.6,2.9,3.2)
z2=c(9.41,5.24,13.84,41.02)

#平均
m1=mean(z1)
m2=mean(z2)

#距離
h=dist(rbind(u1,u2),method="euclidean")
hd=array(0,dim=c(length(z1)+length(z2),length(z1)+length(z2)))
hd[lower.tri(hd)]=c(h);hd=hd+t(hd)


#推定点
h0=c(1.5,3)
h01=sqrt(apply(t(t(u1)-h0)^2,1,sum))
h02=sqrt(apply(t(t(u2)-h0)^2,1,sum))

C=hd
C[1:nrow(u1),1:nrow(u1)]=1.64*exp(-(hd[1:nrow(u1),1:nrow(u1)]/2.91)^2)
C[(nrow(u1)+1):nrow(C),1:nrow(u1)]=2.2*exp(-(hd[(nrow(u1)+1):nrow(C),1:nrow(u1)]/3.2)^2)
C[1:nrow(u1),(nrow(u1)+1):nrow(C)]=2.2*exp(-(hd[1:nrow(u1),(nrow(u1)+1):nrow(C)]/3.2)^2)
C[-c(1:nrow(u1)),-c(1:nrow(u1))]=6.5*exp(-(hd[-c(1:nrow(u1)),-c(1:nrow(u1))]/2.4)^2)

C01=1.64*exp(-(h01/2.91)^2)
C02=2.2*exp(-(h02/3.2)^2)
C0=c(C01,C02)

w=solve(C)%*%C0
z=m1+sum(w*c(z1-m1,z2-m2))
sig=1.64-sum(C0*w)

```




```r


#単純クリギング　p.86
#誤差分散を最小にするレンジ、ナゲット、シルを直接計算する

#各座標
x=c(2,1,4,6,5.5)
y=c(4,3,2.5,2,4)
u=cbind(x,y)

#観測ベクトル
z=c(5.3,4.5,4.6,2.9,3.2)

#距離
h=dist(u,method="euclidean")
hd=array(0,dim=c(length(z),length(z)))
hd[lower.tri(hd)]=c(h);hd=hd+t(hd)

#推定点１
h1=c(1.5,3)
h0=sqrt(apply(t(t(u)-h1)^2,1,sum))

m=mean(z)

#理論バリオグラム
#指数モデル

b=1;c=0.5;a=0.5

ite=10000
eta=0.001

for(i in 1:ite){

rhe=(b+(c-b)*(1-exp(-Inf/a)))-(b+(c-b)*(1-exp(-hd/a)))
rhe0=(b+(c-b)*(1-exp(-Inf/a)))-(b+(c-b)*(1-exp(-h0/a)))
we=solve(rhe)%*%rhe0
ZSKe=m+sum((z-m)*we)
sigSKe=2.5*(1-exp(-Inf/5.53))-sum(we*rhe0)
pthie=c(-1,we)
Ce=rbind(c((b+(c-b)*(1-exp(-Inf/a))),rhe0),cbind(rhe0,rhe))
coste=t(pthie)%*%Ce%*%t(t(pthie))

rhea=-(c-b)*hd*exp(-hd/a)/(a^2)
rhea0=-(c-b)*h0*exp(-h0/a)/(a^2)
Cea=rbind(c(0,rhea0),cbind(rhea0,rhea))
da=t(pthie)%*%(Cea)%*%t(t(pthie))

rheb=-(1-(1-exp(-hd/a)))
rheb0=-(1-(1-exp(-h0/a)))
Ceb=rbind(c(0,rheb0),cbind(rheb0,rheb))
db=t(pthie)%*%(Ceb)%*%t(t(pthie))

rhec=1-(1-exp(-hd/a))
rhec0=1-(1-exp(-h0/a))
Cec=rbind(c(1,rhec0),cbind(rhec0,rhec))
dc=t(pthie)%*%(Cec)%*%t(t(pthie))

a=ifelse(as.numeric(a-eta*da)>0,as.numeric(a-eta*da),0.01)
b=ifelse(as.numeric(b-eta*db)>0,as.numeric(b-eta*db),0.00001)
c=ifelse(as.numeric(c-eta*dc)>0,as.numeric(c-eta*dc),0.00002)

print(as.numeric(coste))

}

```

```r


#単純クリギング　p.86
#誤差分散を最小にするレンジ、ナゲット、シルを計算する

#各座標
x=c(2,1,4,6,5.5)
y=c(4,3,2.5,2,4)
u=cbind(x,y)

#観測ベクトル
z=c(5.3,4.5,4.6,2.9,3.2)

#距離
h=dist(u,method="euclidean")
hd=array(0,dim=c(length(z),length(z)))
hd[lower.tri(hd)]=c(h);hd=hd+t(hd)

#推定点１
h1=c(1.5,3)
h0=sqrt(apply(t(t(u)-h1)^2,1,sum))

m=mean(z)

#理論バリオグラム
#ガウスモデル

b=1;c=0.5;a=0.5

ite=10000
eta=0.0001

for(i in 1:ite){

rhg=c-(b+(c-b)*(1-exp(-(hd/a)^2)))
rhg0=c-(b+(c-b)*(1-exp(-(h0/a)^2)))
wg=solve(rhg)%*%rhg0
ZSKg=m+sum((z-m)*wg)
pthig=c(-1,wg)
Cg=rbind(c(c,rhg0),cbind(rhg0,rhg))
costg=t(pthig)%*%Cg%*%t(t(pthig))

rhea=2*(c-b)*(hd^2)*exp(-(hd/a)^2)/(a^3)
rhea0=2*(c-b)*(h0^2)*exp(-(h0/a)^2)/(a^3)
Cea=rbind(c(0,rhea0),cbind(rhea0,rhea))
da=t(pthie)%*%(Cea)%*%t(t(pthie))

rheb=-(1-(1-exp(-(hd/a)^2)))
rheb0=-(1-(1-exp(-(h0/a)^2)))
Ceb=rbind(c(0,rheb0),cbind(rheb0,rheb))
db=t(pthie)%*%(Ceb)%*%t(t(pthie))

rhec=1-(1-exp(-(hd/a)^2))
rhec0=1-(1-exp(-(h0/a)^2))
Cec=rbind(c(1,rhec0),cbind(rhec0,rhec))
dc=t(pthie)%*%(Cec)%*%t(t(pthie))

a=ifelse(as.numeric(a-eta*da)>0,as.numeric(a-eta*da),0.01)
b=ifelse(as.numeric(b-eta*db)>0,as.numeric(b-eta*db),0.00001)
c=ifelse(as.numeric(c-eta*dc)>0,as.numeric(c-eta*dc),0.00002)

print(as.numeric(costg))

}

```