Edzer and Ashton
Thank you for the tips, i will go back to the drawing board and correct my
mistakes, and sorry everybody for the inconvenience. I used solve to find the
inverse and to get lu. Now i will try to figure out how to use loops to make
the code generic the do as you suggested.
Kabeli
--- On Mon, 10/26/09, Edzer Pebesma <edzer.pebe...@uni-muenster.de> wrote:
From: Edzer Pebesma <edzer.pebe...@uni-muenster.de>
Subject: Re: [R-sig-Geo] Univeral kriging
To: "TANKISO THEJANE" <tankisothej...@yahoo.com>
Cc: R-sig-Geo@stat.math.ethz.ch
Date: Monday, October 26, 2009, 1:34 PM
Tankiso,
Here is some tips:
- make your code such that cut and paste works for us; for me, the ^2
below is vanished (and interpreted as superscript right away, I don't
know why)
- you didn't test your code, as the function inv does not exist
- your code doesn't compute the kriging variance
- the code is very hard to read, and uncommented, so we have to guess
what you want (first order linear trend in the coordinates?)
- use functions if you want to make something generic
- write and use functions such that matrices are allowed, e.g. for g()
you can pass a matrix vector; don't write loops over the individual elements
- If you want to solve Ax=b, write solve(A,b) instead of solve(A) %*% b.
- if you want to find out if your code works, don't post it here but
compare your output with that obtained from packages that do universal
kriging such as fields, geoR or gstat.
--
Edzer
TANKISO THEJANE wrote:
> R geo-helpers
>
> i am new to r and would like you to help me as i am trying to go through
> basic steps. I need to compute a predicted value at point (3,3) and its
> variance given la set of locations and value. i want to make the code generic
> for any size such that i would be able to estimate universal krige at each
> point of grid (0,..,6) (0,..6).
>
> the code
>
> g = function(h) { return(0.1+3*(1-exp(-h/4)))}## variogram of exponential form
> x = rbind(c(1,1),c(2,5),c(4,1),c(5,4))##4 observed locations
> x = rbind(x,c(3,3)) # Row five is the point at which we want to estimate.
> d = sqrt(x[,1]^2+x[,2]^2)##distance of the locations from the origin
> s = c(1,9,5,12)## values at the corresponding 4 locations
> ## pairsv = pair semivariogram
> pairsv = function(i,j) {
> p1 = x[i,]
> p2 = x[j,]
> d = sqrt(sum((p1-p2)^2))
> return(g(d))
> }
> G = matrix(0,4,4)
> for(i in 1:4) for(j in 1:4) G[i,j] = pairsv(i,j)
> X = matrix(0,4,3)
> for(i in 1:4) for(j in 1:3) X[i,j] = d[i]^(j-1)
> ##Asymmetric matrix
> Gu = rbind(cbind(G,X), cbind(t(X),matrix(0,3,3)))
> x
> gu = (1:7)*0
> gu
> for(i in 1:4) gu[i] = pairsv(i,5)
> for(i in 1:3) gu[i+4] = d[5]^(i-1)
> lu = inv(Gu) %*% gu
>
>
>
>
>
>
> [[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, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
http://www.springer.com/978-0-387-78170-9 e.pebe...@wwu.de
__________________________________________________
[[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
