Dear R people,
I have a mysterious problem I have been tearing my hair over for some
time. I have a function that looks as follows:
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)
tempmatrix
}
This returns a matrix of dim pointnum times pointnum with
all entries zero. Now, if I try covexp(topo.point,topo.pair,c(1,1,1),2,3)
I get a matrix of all zeros, which cannot be correct. However, if I do
expbinsumsq(topo.point,topo.pair,1,2,2,3,c(1,1,1)) I don''t get zero
but
some decimal. So at least the (2,3) element of the matrix should have been
non-zero.
Note that my function tempexpbinsumsq merely exists for the purpose of
outer. I could have done
tempmatrix <-
outer(x,y,function(x,y) expbinsumsq(point,pair,x,y,a,b,theta))
but I don''t know if that would be correct usage.
Can someone explain what is going on? I''m at my wits end.
Sincerely, Faheem Mitha.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
At 07:37 PM 5/6/00 -0400, Faheem Mitha wrote:>Dear R people, > >I have a mysterious problem I have been tearing my hair over for some >time. I have a function that looks as follows: > >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) > tempmatrix >}For this to work your function expbinsumsq has to be *vectorized* with respect to x and y. That''s something you have to check. This is the most common reason (in my experience at least) for outer() and other such functions not to perform as expected. (Note, too, that this would have no show of working in S (as opposed to R) because of visibility restrictions. If you were programming for portability this would be an issue.)> >This returns a matrix of dim pointnum times pointnum with >all entries zero. Now, if I try covexp(topo.point,topo.pair,c(1,1,1),2,3) >I get a matrix of all zeros, which cannot be correct. However, if I do >expbinsumsq(topo.point,topo.pair,1,2,2,3,c(1,1,1)) I don''t get zero but >some decimal. So at least the (2,3) element of the matrix should have been >non-zero. > >Note that my function tempexpbinsumsq merely exists for the purpose of >outer. I could have done>tempmatrix <- >outer(x,y,function(x,y) expbinsumsq(point,pair,x,y,a,b,theta))>but I don''t know if that would be correct usage.In R, yes, if your work-horse function expbinsumsq is vectorized. In S you would also need to do something to make the additional arguments available.> >Can someone explain what is going on? I''m at my wits end.Yes, but it would take quite some space here. I''d say get to understand as much as you can about vectorization first, though. I suspect that is your problem.> > Sincerely, Faheem Mitha. > >-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.->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 >_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._> >-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Faheem Mitha <faheem at email.unc.edu> writes:> Note that my function tempexpbinsumsq merely exists for the purpose of > outer. I could have done > tempmatrix <- > outer(x,y,function(x,y) expbinsumsq(point,pair,x,y,a,b,theta)) > but I don''t know if that would be correct usage. > > Can someone explain what is going on? I''m at my wits end.The thing that usually tricks beginners with outer() is that the function has to be vectorised. I.e. if you stick in vectors for x and y, you get a vector result back. Otherwise, you have to vectorise it yourself, e.g. if f takes scalar arguments, f.vect <- function(x,y) sapply(seq(along=x),function(i)f(x[i],y[i])) -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /''_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On 7 May 2000, Peter Dalgaard BSA wrote:> Faheem Mitha <faheem at email.unc.edu> writes: > > > Note that my function tempexpbinsumsq merely exists for the purpose of > > outer. I could have done > > tempmatrix <- > > outer(x,y,function(x,y) expbinsumsq(point,pair,x,y,a,b,theta)) > > but I don''t know if that would be correct usage. > > > > Can someone explain what is going on? I''m at my wits end. > > The thing that usually tricks beginners with outer() is that the > function has to be vectorised. I.e. if you stick in vectors for x and > y, you get a vector result back. Otherwise, you have to vectorise it > yourself, e.g. if f takes scalar arguments, > > f.vect <- function(x,y) sapply(seq(along=x),function(i)f(x[i],y[i]))This is helpful. However I''d like to clarify the meaning of vectorisation in this case. Does a vectorised function with two arguments mean that if X, Y vectors then f(X,Y)= (f(X_i,Y_i)), ie f(X,Y) is the vector with component f(X_i,Y_i)? This is what appears to be the case from the line beginning f.vect above. In particular, this would force x and y to be the same length. If this is not the case then I am puzzled how f.vect would be what I want. A priori, I could have taken this to mean that f must satisfy the following: If X is a vector, and y is a scalar, then f(X,y) = (f(X_i,y)) ie f(X,y) is the vector with components (f(X_i,y) and similarly for f(x,Y). But this is not what you mean, is it? Please excuse any confusion. Trying to debug stuff always wears me out, so I am not at my best right now. Thanks. Faheem. --> O__ ---- Peter Dalgaard Blegdamsvej 3 > c/ /''_ --- Dept. of Biostatistics 2200 Cph. N > (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 > ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907 >-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
(replying to my own message)
Ok, I took a look around, and what I think vectorised means is that if x
and y are two vectors, then the shorter of the vectors is replicated to be
the length of the longer of the two vectors, ie. the elements of the
shorter vector are repeated as necessary (thus the two vectors become the
same length) and then a vector of function values is returned that is the
value of the function elementwise on these two vectors. (Thanks to S
Burn''s `S Poetry'' for a clear explanation of this on pg 4). I
am still not
clear how I would vectorise my function though, assuming that it only
accepts scalar values.
Cam some kind person in any case confirm that this is what was meant?
f.vect <- function(x,y) sapply(seq(along=x),function(i)f(x[i],y[i]))
seems to be to only work correctly if x and y are the same length.
Sincerely, Faheem Mitha.
On Sat, 6 May 2000, Faheem Mitha wrote:
>
>
> On 7 May 2000, Peter Dalgaard BSA wrote:
>
> > Faheem Mitha <faheem at email.unc.edu> writes:
> >
> > > Note that my function tempexpbinsumsq merely exists for the
purpose of
> > > outer. I could have done
> > > tempmatrix <-
> > > outer(x,y,function(x,y) expbinsumsq(point,pair,x,y,a,b,theta))
> > > but I don''t know if that would be correct usage.
> > >
> > > Can someone explain what is going on? I''m at my wits
end.
> >
> > The thing that usually tricks beginners with outer() is that the
> > function has to be vectorised. I.e. if you stick in vectors for x and
> > y, you get a vector result back. Otherwise, you have to vectorise it
> > yourself, e.g. if f takes scalar arguments,
> >
> > f.vect <- function(x,y)
sapply(seq(along=x),function(i)f(x[i],y[i]))
>
> This is helpful. However I''d like to clarify the meaning of
vectorisation
> in this case.
>
> Does a vectorised function with two arguments mean that if X, Y vectors
> then f(X,Y)= (f(X_i,Y_i)), ie f(X,Y) is the vector with component
> f(X_i,Y_i)?
>
> This is what appears to be the case from the line beginning f.vect above.
> In particular, this would force x and y to be the same length. If this is
> not the case then I am puzzled how f.vect would be what I want.
>
> A priori, I could have taken this to mean that f must satisfy the
> following:
>
> If X is a vector, and y is a scalar, then f(X,y) = (f(X_i,y)) ie f(X,y) is
> the vector with components (f(X_i,y) and similarly for f(x,Y). But this is
> not what you mean, is it?
>
> Please excuse any confusion. Trying to debug stuff always wears me out, so
> I am not at my best right now.
>
> Thanks. Faheem.
>
> --
> > O__ ---- Peter Dalgaard Blegdamsvej 3
> > c/ /''_ --- Dept. of Biostatistics 2200 Cph. N
> > (*) \(*) -- University of Copenhagen Denmark Ph: (+45)
35327918
> > ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45)
35327907
> >
>
>
>
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
At 10:57 PM 5/6/00 -0400, Faheem Mitha wrote:>(replying to my own message) > >Ok, I took a look around, and what I think vectorised means is that if x >and y are two vectors, then the shorter of the vectors is replicated to be >the length of the longer of the two vectors, ie. the elements of the >shorter vector are repeated as necessary (thus the two vectors become the >same length) and then a vector of function values is returned that is the >value of the function elementwise on these two vectors. (Thanks to S >Burn''s `S Poetry'' for a clear explanation of this on pg 4).It is a clear, concise and simple explanation, but sadly incomplete and for your purposes rather misleading... Darn!>I am still not clear how I would vectorise my function though, assuming >that it only accepts scalar values.How you do that is up to you. There is no magic formula. Sometimes efficiency requires you to go right down to the C code (or extension) level and write a fast loop.> >Cam some kind person in any case confirm that this is what was meant? > >f.vect <- function(x,y) sapply(seq(along=x),function(i)f(x[i],y[i])) > >seems to be to only work correctly if x and y are the same length.Call me unkind, but this is precisely the kind of discussion I was seeking to avoid in public as it gets so messy. It takes a dialogue. "Vectorization" is a rather slippery concept that means different things in different contexts. It may involve the recycling rule or it may not, depending. Take an example. If I write f <- function(x, y) x + y this function will work with outer(), but it depends on a rather different vectorization property than just the recycling rule. In this case it works because if either x or y are vectors (but not both), the result is a vector of that length. Note that if I use f(1:3, 2) or f(2, 1:3) the recycling rule is used, but it is not used in, say, outer(1:3, 0:5, f) (which you should check). So your function to be supplied to outer() should> Sincerely, Faheem Mitha. > >On Sat, 6 May 2000, Faheem Mitha wrote: > >> >> >> On 7 May 2000, Peter Dalgaard BSA wrote: >> >> > Faheem Mitha <faheem at email.unc.edu> writes: >> > >> > > Note that my function tempexpbinsumsq merely exists for the purpose of >> > > outer. I could have done >> > > tempmatrix <- >> > > outer(x,y,function(x,y) expbinsumsq(point,pair,x,y,a,b,theta)) >> > > but I don''t know if that would be correct usage. >> > > >> > > Can someone explain what is going on? I''m at my wits end. >> > >> > The thing that usually tricks beginners with outer() is that the >> > function has to be vectorised. I.e. if you stick in vectors for x and >> > y, you get a vector result back. Otherwise, you have to vectorise it >> > yourself, e.g. if f takes scalar arguments, >> > >> > f.vect <- function(x,y) sapply(seq(along=x),function(i)f(x[i],y[i])) >> >> This is helpful. However I''d like to clarify the meaning of vectorisation >> in this case. >> >> Does a vectorised function with two arguments mean that if X, Y vectors >> then f(X,Y)= (f(X_i,Y_i)), ie f(X,Y) is the vector with component >> f(X_i,Y_i)? >> >> This is what appears to be the case from the line beginning f.vect above. >> In particular, this would force x and y to be the same length. If this is >> not the case then I am puzzled how f.vect would be what I want. >> >> A priori, I could have taken this to mean that f must satisfy the >> following: >> >> If X is a vector, and y is a scalar, then f(X,y) = (f(X_i,y)) ie f(X,y) is >> the vector with components (f(X_i,y) and similarly for f(x,Y). But this is >> not what you mean, is it? >> >> Please excuse any confusion. Trying to debug stuff always wears me out, so >> I am not at my best right now. >> >> Thanks. Faheem. >> >> -- >> > O__ ---- Peter Dalgaard Blegdamsvej 3 >> > c/ /''_ --- Dept. of Biostatistics 2200 Cph. N >> > (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 >> > ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907 >> > >> >> >> > >-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.->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 >_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._> >-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._