I regret to say I have to send my last message again as it
accidently became polluted by MicroSoft junk. The most
pernicious garbling was to replace the innocent <- by the
egregious <<- globally. Now you can see why I MUST send a
corrected version!
Here is what I really wanted to say below. My apologies for the
disruption.
Bill Venables.
-----------------------------------------------------------------
At 17:33 15/09/00 -0400, Joseph S. Verducci wrote:
>
> There seems to be a minor problem with the first solution
> and the second solution is limited by the stacking of
> the recursion. For standard settings, the limit is 15 columns:
>
> > unix.time(Subsets(15,2))
> [1] 0.48 0.01 0.00 0.00 0.00
>
> > unix.time(Subsets(16,2))
> Error: Expressions nested beyond limit (256)
> Timing stopped at: 2.41999999999996 0.00999999999999979 3 0 0
No, I would not regard those functions as being at all
appropriate here. They were written for quite different
purposes.
>
> Combining modifications of both solutions seems to work
> reasonable well:
>
> Modifying Gardar's function --
>
> Pairs2 <- function(n){
> #
> # returns an {n(n+1)/2} x 2 matrix of paired integers
> #
> a <- matrix(0,n,n)
> nr <- row(a)
> nc <- col(a)
> ut <- !lower.tri(a)
> cbind(nc[ut],nr[ut])}
>
> Then using Steve's idea
>
> PamG <- function(x){
> #
> # returns an {nc(nc+1)/2} x nr} matrix whose columns
> # give the required products
> #
> apply(Pairs2(ncol(x)),1,function(y) {x[,y[1]] * x[,y[2]]})}
Despite appearances this is really a loop solution. I suspect
your best solution is going to be an even more explicit loop,
though. Here's my suggestion:
Quad <- function(X) {
nc <- ncol(X)
X2 <- matrix(0, nrow(X), (nc * (nc + 1))/2)
m <- 1
for(j in 1:nc) {
X2[, m:(m + nc - j)] <- X[, j] * X[, j:nc]
m <- m + nc - j + 1
}
X2
}
Call as
X2 <- Quad(X)
of course.
This works if X is a data frame or a matrix, but the value, X2,
is always a matrix. If speed were a consideration I would
suggest you work with a simple, unadorned matrices X and X2,
though, (i.e. no names or dimnames attributes).
The most important line here is the one setting up the initial
X2, namely
X2 <- matrix(0, nrow(X), (nc * (nc + 1))/2)
It sets aside enough storage to hold the final answer, so you
need only make one trip to the allocator. The loop itself is
then a very tiny overhead.
Bill Venables.
--
Bill Venables, Statistician Tel. +61 7 3826 7251
CSIRO Marine Laboratories, Fax. +61 7 3826 7304
Cleveland, Qld, 4163 Email: Bill.Venables@cmis.csiro.au
AUSTRALIA http://www.cmis.csiro.au/bill.venables/
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
|