More than 1 year has passed since last update.

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

RStudio

Last updated at 2023-08-04Posted at 2021-10-01

#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))
}


#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))


#通常クリギング

#各座標
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)



#平均値クリギング

#各座標
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


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

#各座標
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))


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

#各座標
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))


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

#各座標
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)



#単純クリギング　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))

}



#単純クリギング　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))

}

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up

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

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