[R] problems with outer

[Faheem Mitha]

Dear R people,

I'm taking the liberty of sending out my code. I am using the data set
topo from the V&R package spatial and the data structures point and pair
from the package sgeostat.

It may look a bit long, but most of it is probably not relevant.

I feel reasonably sure that the expbinsumsq function (towards the bottom)
is causing the problem with the vectorisation and outer. I am using an
ifelse condition at the beginning, and that is probably causing the
trouble, since the values I get are precisely what I would get if the
condition was never true. I should mention that the a, b's mentioned are
factors, taking values 1,2...10.

Perhaps someone can suggest a painless way to sort it out.

Thanks.                                               Faheem. 

**********************************************************************

data(topo)
# defining the point object (see package `sgeostat').
topo.point <- point(topo)
#defining the pair object (see the package `sgeostat').
topo.pair <- pair(topo.point)
#Estimating values of the exponential variogram
topo.est.var <- est.variogram(topo.point,topo.pair,'z')


#Making a vector of n(h_j)'s. (The number of pairs in each bin)
binnum <- function(pair)
{ numofbins <- length(pair$bins)
  binnumval <- rep(0,numofbins)
  for(i in 1:numofbins)
{  binnumval[i] <- sum(ifelse(as.vector(pair$lags)==i,1,0))}
   return(binnumval)  
}

# distance between the points (a1,a2) and (b1,b2).
metric <- function(a1,a2,b1,b2)
{ sqrt((a1 - b1)^2 + (a2 - b2)^2) }

# the exponential variogram function
expvar <- function(t,theta)
{
  theta[1] + theta[2]*(1 - exp(-t/theta[3]))
}

expsumsq <- function(x1,y1,x2,y2,x3,y3,x4,y4,theta)
{ ( expvar(metric(x1,y1,x3,y3),theta) + expvar(metric(x2,y2,x4,y4),theta)
  - expvar(metric(x1,y1,x4,y4),theta) - expvar(metric(x2,y2,x3,y3),theta) )^2
}

expbinsumsq <- function(point,pair,i,j,a,b,theta)
{
  ifelse((as.integer(pair$lags[i]) == a)&&(as.integer(pair$lags[j]) == b),
            expsumsq(point$x[pair$from[i]],point$y[pair$from[i]],
                   point$x[pair$to[i]],point$y[pair$to[i]],
                   point$x[pair$from[j]],point$y[pair$from[j]],
                   point$x[pair$to[j]],point$y[pair$to[j]],
                   theta),0)
}

covexp <- function(point,pair,theta,a,b)
{ pointnum <- length(pair$lags)
  
  x <- seq(pointnum)
  y <- x
  tempexpbinsumsq <- function(x,y) {expbinsumsq(point,pair,x,y,a,b,theta)}
  tempmatrix <- outer(x,y,tempexpbinsumsq)
  tempsum <- sum(tempmatrix)
  return( 2*tempsum/( binnum(pair)[a]*binnum(pair)[b] ) )
}

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._