[R] Re: geoR size limit problem

[Ole F. Christensen]

Dear Miha

I can reproduce your problem. The cause of the problem is not the size 
of the data (2500), but a combination of the number data locations and 
prediction location [could be that the covariance matrix of the data and 
the prediction variables is the computational bottleneck].
A workaround by splitting the set of prediction locations into smaller 
blocks is given below.

As I wrote in my previous e-mail, I would consider trying the kriging 
functions in the geostatistical packages.


## Example

library(geoR)
g.data <-as.geodata(cbind(rnorm(2500),rnorm(2500), rnorm(2500)))

expvar<-variog(g.data, uvec=seq(0,1000,25), option="bin")
maternfit<-variofit(expvar, ini=c(1600,300), cov.model="matern", 
fix.nug=TRUE, nug=25, kappa=1.5, max.dist=800, weights="npairs")


p.grid<-expand.grid((1:900)/900,(1:650)/650)

## this doesn't work
krige.matern<-krige.conv(g.data, loc=p.grid, 
krige=krige.control(obj.m=maternfit, type="OK"))


## workaround (be patient, it is slow):
nm.grid <- nrow(p.grid)
result <- list(predict=rep(0,nm.grid), krige.var=rep(0,nm.grid))
claas(result) <- "kriging"
##for(ii in 1:585){
for(ii in 1:5){
temp <- krige.matern<-krige.conv(g.data, 
loc=p.grid[(ii-1)*1000+(1:1000),], krige=krige.control(obj.m=maternfit, 
type="OK"))[c("predict", "krige.var")]
result$predict[(ii-1)*1000+(1:1000)] <- temp$predict
result$krige.var[(ii-1)*1000+(1:1000)] <- temp$krige.var
print(ii)
}


Best Regards Ole



Miha STAUT wrote:

>> Dear Miha (cc: Paulo Ribeiro [the developer of geoR, who is not 
>> subscribed to R-help])
>>
>> Some reproducible code would help here.
>
>
> Unfortunately I do not have an internet connect on my computer, that 
> is why I did every thing very descriptive. Anyway I will try:
>
> str(df)
>   [#approx 2500]
> $ x ...
> $ y ...
> $ z ...
> library(geoR)
> g<-as.geodata(df)
> expvar<-variog(g, uvec=seq(0,1000,25), option="bin")
> matern<-variofit(expvar, ini=c(1600,300), cov.model="matern", 
> fix.nug=T, nug=25, kappa=1.5, max.dist=800, weights="npairs")
>
> library(GRASS)
> G<-gmeta()
> grid<-expand.matrix(G$xseq,G$yseq)
> length(G$xseq)
> 900 #approx
> length(G$yseq)
> 650 #approx
> krige.matern<-krige.conv(g, loc=grid, 
> krige=krige.control(obj.m=matern, type="OK"))
>
> #5 min processing
> can not allocate vector of 1500000000 #approx
>
> Thanks, Miha Staut
>
>>
>> Guessing :
>>
>> * geoR does not implement kriging with local neighbourhoods, but 
>> instead conditions on all data. Therefore having 2500 data points 
>> would imply having a covariance matrix of size 2500*2500 approx 
>> 6million . Not sure if this is too much, but I do not think so.
>
>
> _________________________________________________________________

> http://join.msn.com/?page=features/virus
>
>

-- 
Ole F. Christensen
Center for Bioinformatik
Datalogisk Institut
Aarhus Universitet
Ny Munkegade, Bygning 540
8000 Aarhus C
Denmark