Go to Qiita Advent Calendar 2024 Top
LoginSignup
2
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

@kozakai-ryouta(亮太 小酒井)

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

Last updated at Posted 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))

}
2
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?