R 中 gstat 包中 krige 花费的时间
Time taken to krige in gstat package in R
以下 R 程序使用 gstat 包中的 walker Lake 数据使用 470 个数据点创建一个插值表面。
source("D:/kriging/allfunctions.r") # Reads in all functions.
source("D:/kriging/panel.gamma0.r") # Reads in panel function for xyplot.
library(lattice) # Needed for "xyplot" function.
library(geoR) # Needed for "polygrid" function.
library(akima)
library(gstat);
library(sp);
walk470 <- read.table("D:/kriging/walk470.txt",header=T)
attach(walk470)
coordinates(walk470) = ~x+y
walk.var1 <- variogram(v ~ x+y,data=walk470,width=10) #the width has to be tuned resulting different point pairs
plot(walk.var1,xlab="Distance",ylab="Semivariance",main="Variogram for V, Lag Spacing = 5")
model1.out <- fit.variogram(walk.var1,vgm(70000,"Sph",40,20000))
plot(walk.var1, model=model1.out,xlab="Distance",ylab="Semivariance",main="Variogram for V, Lag Spacing = 10")
poly <- chull(coordinates(walk470))
plot(coordinates(walk470),type="n",xlab="X",ylab="Y",cex.lab=1.6,main="Plot of Sample and Prediction Sites",cex.axis=1.5,cex.main=1.6)
lines(coordinates(walk470)[poly,])
poly.in <- polygrid(seq(2.5,247.5,5),seq(2.5,297.5,5),coordinates(walk470)[poly,])
points(poly.in)
points(coordinates(walk470),pch=16)
coordinates(poly.in) <- ~ x+y
krige.out <- krige(v ~ 1, walk470,poly.in, model=model1.out)
print(krige.out)
本程序对2688个点的每个点计算如下
(470x470) matrix inversion
(470x470) and (470x1) matrix multiplication
gstat 包是否使用了一些智能的计算方法。我从以前的 Whosebug 查询中知道它使用 cholesky 分解进行矩阵求逆。一台机器算这么快算正常吗
它使用LDL'分解，类似于Choleski。当您使用全局克里金法时，协方差矩阵只需要分解一次；然后，对于每个预测点，求解一个系统，即O(n)。没有 470x470 矩阵被反转，也没有通过乘法得到的解决方案。逆是符号设备，但尽可能避免作为计算策略。例如，在 R 中，比较 solve(A,b) 与 solve(A) %*% b.
的运行时间
使用source，卢克！
以下 R 程序使用 gstat 包中的 walker Lake 数据使用 470 个数据点创建一个插值表面。
source("D:/kriging/allfunctions.r") # Reads in all functions.
source("D:/kriging/panel.gamma0.r") # Reads in panel function for xyplot.
library(lattice) # Needed for "xyplot" function.
library(geoR) # Needed for "polygrid" function.
library(akima)
library(gstat);
library(sp);
walk470 <- read.table("D:/kriging/walk470.txt",header=T)
attach(walk470)
coordinates(walk470) = ~x+y
walk.var1 <- variogram(v ~ x+y,data=walk470,width=10) #the width has to be tuned resulting different point pairs
plot(walk.var1,xlab="Distance",ylab="Semivariance",main="Variogram for V, Lag Spacing = 5")
model1.out <- fit.variogram(walk.var1,vgm(70000,"Sph",40,20000))
plot(walk.var1, model=model1.out,xlab="Distance",ylab="Semivariance",main="Variogram for V, Lag Spacing = 10")
poly <- chull(coordinates(walk470))
plot(coordinates(walk470),type="n",xlab="X",ylab="Y",cex.lab=1.6,main="Plot of Sample and Prediction Sites",cex.axis=1.5,cex.main=1.6)
lines(coordinates(walk470)[poly,])
poly.in <- polygrid(seq(2.5,247.5,5),seq(2.5,297.5,5),coordinates(walk470)[poly,])
points(poly.in)
points(coordinates(walk470),pch=16)
coordinates(poly.in) <- ~ x+y
krige.out <- krige(v ~ 1, walk470,poly.in, model=model1.out)
print(krige.out)
本程序对2688个点的每个点计算如下
(470x470) matrix inversion
(470x470) and (470x1) matrix multiplication
gstat 包是否使用了一些智能的计算方法。我从以前的 Whosebug 查询中知道它使用 cholesky 分解进行矩阵求逆。一台机器算这么快算正常吗
它使用LDL'分解，类似于Choleski。当您使用全局克里金法时，协方差矩阵只需要分解一次；然后，对于每个预测点，求解一个系统，即O(n)。没有 470x470 矩阵被反转，也没有通过乘法得到的解决方案。逆是符号设备，但尽可能避免作为计算策略。例如，在 R 中，比较 solve(A,b) 与 solve(A) %*% b.
使用source，卢克！