Re: trying to use octave-gpr

[Jaroslav Hajek]

On Sun, Jul 18, 2010 at 11:18 AM, Oz Nahum <address@hidden> wrote:
> Hi Everyone,
>
> I've recently discovered octave-gpr package, and I would like to use it.
> Until now I have been using R gstat for doing that kind of tasks - namely,
> kriging for predicting spatial data.
>

Beware that OctGPR was written from the viewpoint of an amateur in
statistical prediction (me), who knows very little about
geostatistics. I will try to help you using the hopeless minimum of
geostatistical terminology I managed to imbibe about the problem.

> According to the terminology I used to work with - I need to create a
> variogram for a data set before I can create a krigging prediction.
> I figured out that creating the variogram is done by the function gpr_train,
> and the actually predicted surface is done by gpr_predict .
>

Not quite. gpr_train won't help you to create an empirical variogram,
nor will it attempt to construct one (and fit a model to it).
Rather, OctGPR fits the parameters of a variogram model by directly
maximizing the marginal likelihood of the observation data.
A more or less fixed variogram model is used in the form
gamma(h) = var * f(norm(theta.*h)) + nu * all(h == 0)

where f is one of the decay functions described in gpr_train's
docstring (default is gaussian). In other words, the random spatial
process is considered stationary but *not* isotropic, with noise
(nugget effect) allowed but decaying to zero (sill = 0). This spatial
process is assumed to be superimposed onto a constant or linear trend
surface.

I'm attaching a paper I wrote about OctGPR some time ago, perhaps it
can clarify more.

> Never the less I failed to use the functions, because I find the
> documentation a bit confusing.
>
> My data is loaded from a file using the following structure:
> "x", "y", "obseravation" columns.
>
> next I want to create a "variogram" using
> octave:30> g=gpr_train(xs,ys,t)
> warning: gpr_train: some elements in list of return values are undefined
>

I assume you want to do (the simplest setup possible, tweaking may
yield better results):

model = gpr ([xs; ys], t, theta, {});

where theta are initial estimates for the spatial scales. In absence
of better guesses, I typically use something like
theta = [1./std(xs), 1./std(ys)];

further parameters may specify the initial noise level, linear trend
(default is constant, i.e. "ordinary kriging") and various parameters
for maximum likelihood optimization in a cell array. If you omit the
cell array, no optimization will take place and the function will just
return the model and its negative log likelihood. This can be used to
employ your own estimation techniques (e.g. an empirical variogram).

> I don't really understand what the values of t should be.
>
> For doing a krigging in R I would do the following:
>
> #define coordinates
> coordinates(sic.val) = ~x+y
> #make variogram
> v<-variogram(dataset$dayx~1,dataset["dayx"])
> #find best fit with 120,200,10 as initial guesses
> fit.Gau.v<-fit.variogram(v, vgm(120, "Gau",200,10))
> #create kriging predictor
> xv.Gau <- krige.cv(dayx~1, sic.val["dayx"], fit.Gau.v, nmax = 40, nfold=5)
>
> Can someone please give me an example how to do the same thing with the
> gpr_package ?

Unfortunately I have no idea what the above commands do; so this is a
dead end :)

> I was trying to find the demos here
> but the link seems to be  dead
> :http://artax.karlin.mff.cuni.cz/~hajej2am/octgpr.php

It works normally for me. The server is definitely running.

You can also try running
demo_octgpr (1)

and then look at the source to see how things are done.

best regards

-- 
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

m6-hajek.pdf
Description: Adobe PDF document