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

#p.74 理論バリオグラム(最小二乗法でのあてはめ)#実際に当てはめのデータを作成する際には近しい区間で区切って#値を平均化したものを使用するdata(EuStockMarkets)u=EuStockMarkets[1:200,1:4]hs=array(0,dim=c(nrow(u),nrow(u)))for(iin1:nrow(hs)){hs[i,]=apply((t(u)-c(u[i,]))^2,2,sum)}r=2+(3-2)*(1-exp(-hs/100000))r=r*100Model=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 Modelpredicts=ifelse(h<a,b+(c-b)*(3*h/(2*a)-0.5*(h/a)^3),c)}if(model==2){#Exponential Modelpredicts=b+(c-b)*(1-exp(-h/a))}if(model==3){#Gaussian Modelpredicts=b+(c-b)*(1-exp(-(h/a)^2))}if(model==4){#Hole-Effect Cyclic Modelpredicts=b+(c-b)*(1-cos(pi*h/a))}return(predicts)}f=function(x){vec=c(c(hs),x)predicts=Model(vec)return(predicts)}model=1Y=c(r)h=0.001lam=0.01ite=10^4sita=c(mean(hs[hs>0]),rep(0.01,2))eta=10^(-2)for(lin1:ite){sita_pre=sitaZ=array(0,dim=c(length(Y),length(sita)))vec=sitafor(jin1:ncol(Z)){vec_sub=vec;vec_sub[j]=vec_sub[j]+hZ_sub=(f(vec_sub)-f(vec))/hZ[,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.79r=dist(z,method="euclidean")rd=array(0,dim=c(length(z),length(z)))rd[lower.tri(rd)]=c(r);rd=rd+t(rd)/2plot(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]=0rhk0=1.84-1.84*(3*h0/(2*7.69)-(1/2)*(h0/7.69)^3)rhk0[h0>7.69]=0wk=solve(rhk)%*%rhk0ZSKk=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)%*%rhg0ZSKg=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)%*%rhe0ZSKe=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.79r=dist(z,method="euclidean")rd=array(0,dim=c(length(z),length(z)))rd[lower.tri(rd)]=c(r);rd=rd+t(rd)/2plot(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.84rhk=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.84wk=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.79r=dist(z,method="euclidean")rd=array(0,dim=c(length(z),length(z)))rd[lower.tri(rd)]=c(r);rd=rd+t(rd)/2plot(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.84rhk=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)%*%zf=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.79r=dist(z,method="euclidean")rd=array(0,dim=c(length(z),length(z)))rd[lower.tri(rd)]=c(r);rd=rd+t(rd)/2plot(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]=0rhk=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]=0wk=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)*bS=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]=0rhk=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]=0wk=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=hdC[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)%*%C0z=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.5ite=10000eta=0.001for(iin1: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)%*%rhe0ZSKe=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.5ite=10000eta=0.0001for(iin1:ite){rhg=c-(b+(c-b)*(1-exp(-(hd/a)^2)))rhg0=c-(b+(c-b)*(1-exp(-(h0/a)^2)))wg=solve(rhg)%*%rhg0ZSKg=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))}