Edzer, Thanks for identifying the problem. The problem showed how big the impact of the choice of nmax value on the estimation is.
In both krige and gstat, the default value for nmax is Inf. I am doing a simulation experiment now to assess the performance of several spatial interpolators using a few hundreds datasets. It is easy to do it by using the default value. But I am wondering what the best guess for nmax is? I assumed that the default value Inf for nmax would take the maximum value of the number of samples available. After a few trials, I found my assumption was wrong as evidenced below. Obviously, it took a value between 20 and 30. > system.time(x <- krige(log(zinc)~1, meuse, meuse.grid, model = m)) [using ordinary kriging] user system elapsed 0.32 0.00 0.31 > system.time(x <- krige(log(zinc)~1, meuse, meuse.grid, model = m, nmax=20)) [using ordinary kriging] user system elapsed 0.25 0.02 0.27 > system.time(x <- krige(log(zinc)~1, meuse, meuse.grid, model = m, nmax=30)) [using ordinary kriging] user system elapsed 0.41 0.02 0.43 > system.time(x <- krige(log(zinc)~1, meuse, meuse.grid, model = m, nmax=50)) [using ordinary kriging] user system elapsed 0.89 0.00 0.89 Please clarify this. Thanks, Jin -----Original Message----- From: Edzer Pebesma [mailto:[EMAIL PROTECTED] Sent: Tuesday, 8 July 2008 5:35 To: Li Jin Cc: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] why are idw and gstat different? [SEC=UNCLASSIFIED] Jin, you specified nmax = 7 in the second call, but not in the first. Comparing idp values is easiest when other things remain equal. -- Edzer [EMAIL PROTECTED] wrote: > Dear there, > > I tried to compare idw and gstat in library(gstat). I specified idp as 0.5 in > both functions, I expected identical results, but what I got are different as > shown for the first five samples. Please help. Thanks. > > > > >> data(meuse) >> > > >> coordinates(meuse) = ~x+y >> > > >> data(meuse.grid) >> > > >> gridded(meuse.grid) = ~x+y >> > > >> x <- idw(zinc~1, meuse, meuse.grid, idp=0.5) >> > > [inverse distance weighted interpolation] > > >> as.data.frame(x)[1:5,] >> > > var1.pred var1.var x y > > 1 482.8753 NA 181180 333740 > > 2 487.0065 NA 181140 333700 > > 3 483.9747 NA 181180 333700 > > 4 481.2115 NA 181220 333700 > > 5 494.8720 NA 181100 333660 > > > > >> data(meuse) >> > > >> data(meuse.grid) >> > > >> meuse.gstat <- gstat(id = "zinc", formula = zinc ~ 1, locations = ~ x + y, >> > > + data = meuse, nmax = 7, set = list(idp = .5)) > > >> meuse.gstat >> > > data: > > zinc : formula = zinc`~`1 ; data dim = 155 x 12 nmax = 7 > > set idp = 0.5; > > ~x + y > > >> z <- predict(meuse.gstat, meuse.grid) >> > > [inverse distance weighted interpolation] > > >> z[1:5,] >> > > x y zinc.pred zinc.var > > 1 181180 333740 626.3628 NA > > 2 181140 333700 645.9319 NA > > 3 181180 333700 629.7041 NA > > 4 181220 333700 615.1368 NA > > 5 181100 333660 682.3401 NA > > > > Cheers, > > > > Jin > > -------------------------------------------- > > Jin Li, PhD > > Spatial Modeller/ > > Computational Statistician > > Marine & Coastal Environment > > Geoscience Australia > > > > Ph: 61 (02) 6249 9899 > > Fax: 61 (02) 6249 9956 > > email: [EMAIL PROTECTED] <mailto:[EMAIL PROTECTED]> > > -------------------------------------------- > > > > > [[alternative HTML version deleted]] > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@stat.math.ethz.ch > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics (IfGI), University of Münster, Weseler Straße 253, 48151 Münster. Phone: +49 251 8333081 Fax: +49 251 8339763 email: [EMAIL PROTECTED] http://ifgi.uni-muenster.de/ _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
