Hi Zev,
There's no need to brute force, as optimize is there to help you -- my
guess is that the function is convex. The following takes a while:
> f = function(idp, formula, data,...)
sum(krige.cv(formula,data,set=list(debug=0,idp=idp),...)$residual**2)
> optimize(f, interval=c(0.01,4), formula=log(zinc)~1, data=meuse)
$minimum
[1] 3.532051
$objective
[1] 32.09126
but that is probably due to the (my?) hopelessly inefficient (though
flexible) setup of krige.cv. Speeding up can be done by allowing a
larger tolerance, or passing e.g. nfold=5 to optimize(). This nfold will
also result in somewhat random output, as it randomly folds the data in
5 partitions.
--
Edzer
Zev Ross wrote:
Hi All,
Thanks to Paul and Alessandro for their suggestions. Paul's code
(brute force) worked well for me and the results match up well with
ArcGIS. I'm not using a large dataset so the speed isn't an issue but
with a larger dataset it would be. In ArcGIS the optimization is
instantaneous so clearly the software is doing something different
than testing out all possible values.
I used Paul's code with no substantive changes then it's
straightforward to use:
plot(idw_pow, cv_vals)
idw_pow[which.min(cv_vals)]
And pull out the min.
Thanks for the advice! Zev
Alessandro wrote:
Hi
Normally I use the R+SAGA to calculate the IDW and create a raster, with
this follow code. I change the radius input with a loop formula to create
several raster and after check the best result (I am studying a oak forest
with low density)
radii.list <- as.list(c(5, 10, 15, 20, 25, 30))
for(i in 1:length(radii.list)){
rsaga.geoprocessor(lib="grid_gridding", module=0,
param=list(GRID=set.file.extension(paste("DEM_iw",radii.list[[i]],sep=""),".
sgrd"),
SHAPES="ground.shp", FIELD=0, TARGET= 0, POWER=1.0,
RADIUS=radii.list[[i]], NPOINTS=20, USER_CELL_SIZE=dem.pixelsize,
[EMAIL PROTECTED],1], [EMAIL PROTECTED],2],
[EMAIL PROTECTED],1], [EMAIL PROTECTED],2]))
}
After I have 6 raster (DEM_idw_5.sgrd, DEM_idw_10.sgrd, DEM_idw_15.sgrd,
DEM_idw_20.sgrd, etc. etc.). I check hand-made the best IDW. I don't know is
It possible to improve this formula with RMSE in gstat and calculate the
best power.
Ale
-----Messaggio originale-----
Da: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] Per conto di Paul Hiemstra
Inviato: martedì 2 dicembre 2008 12.07
A: Zev Ross
Cc: r-sig-geo@stat.math.ethz.ch
Oggetto: Re: [R-sig-Geo] GSTAT: Optimize power value for IDW
Zev Ross schreef:
Hi All,
ArcGIS has a nice little button that calculates the optimal power
value to use for inverse distance weighting based on cross-validation
and RMSPE. Just wondering if anyone had written something similar in R
-- I'm using GSTAT and I'd like to avoid back and forth with ArcGIS
(and obviously I'd like to avoid writing it myself as well!).
Zev
Hi,
I don't have any code lying around, but a brute force optimization
approach should be quite easy. Also because the range of idw-powers is
relatively small. The speed would ofcourse depend on the size of the
dataset. In code it would look something like (actually works :)):
library(gstat)
data(meuse)
coordinates(meuse) = ~x+y
# make function to do the cv
do_cv = function(idp) {
meuse_idw = gstat(id = "zn", formula = zinc~1, data = meuse, set =
list(idp = idp))
out = gstat.cv(meuse_idw)
return(sqrt(mean(out$residual^2))) # RMSE as optimization criterium
}
idw_pow = seq(0.1,5, by = 0.2) # the idwpower values that will be checked
cv_vals = sapply(idw_pow, do_cv) # calculate the rmse
# List of outcomes
print(data.frame(idp = idw_pow, cv_rmse = cv_vals))
After this you select the idw value with the smallest RMSE.
cheers and hth,
Paul
--
Zev Ross
ZevRoss Spatial Analysis
120 N Aurora, Suite 3A
Ithaca, NY 14850
607-277-0004 (phone)
866-877-3690 (fax, toll-free)
[EMAIL PROTECTED]
------------------------------------------------------------------------
_______________________________________________
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, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
http://www.springer.com/978-0-387-78170-9 [EMAIL PROTECTED]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
