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
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]]})}
___________________________________________________
On Fri, 15 Sep 2000, Pamela Goodman wrote:
>
>
>
> Thanks to Gardar Johannesson and Steve McKinney^2 for their suggestions. I
> have not tested them for my needs yet, but, they will likely do very well
> for me.
>
>
> ### from Gardar Johannesson
>
> foo <- function(x) {
> xr <- row(x)
> xc <- col(x)
> ind <- xr <= xc
> x[,xr[ind]] * x[,xc[ind]]
> }
>
> x is the matrix. the function returns a matrix with each column storing
> the product of some two columns in the matrix x.
>
>
> ### from Steve McKinney [1]
> This function by Bill Venables might
> do the trick, try it out.
> You can just add the indices such as
> 1 1
> 2 2
> 3 3
> etc.
> to the list produced by subsets().
>
>
> subsets <- function(n, r, v = 1:n)
> if(r <= 0) vector(mode(v), 0) else
> if(r >= n) v[1:n] else {
> rbind(cbind(v[1], Recall(n-1, r-1, v[-1])),
> Recall(n-1, r, v[-1]))
> }
>
>
> > subsets(6,2)
> [,1] [,2]
> [1,] 1 2
> [2,] 1 3
> [3,] 1 4
> [4,] 1 5
> [5,] 1 6
> [6,] 2 3
> [7,] 2 4
> [8,] 2 5
> [9,] 2 6
> [10,] 3 4
> [11,] 3 5
> [12,] 3 6
> [13,] 4 5
> [14,] 4 6
> [15,] 5 6
> >
>
> To avoid extra computations, you could run subsets
> for some large value of n, save the results, and
> just access the saved results to get the index
> info you need.
>
>
> ### from Steve McKinney [2]
>
> Here's another version of Bill Venable's
> code for subset calculations that saves
> intermediate results for speedier processing.
>
>
> Subsets <- function(n, r, v = 1:n, frame = 1)
> if(r <= 0) vector(mode(v), 0) else
> if(r >= n) v[1:n] else {
> if(r > 1) {
> i1 <- paste("#", n - 1, r - 1)
> i2 <- paste("#", n - 1, r)
> if(!exists(i1))
> assign(i1,
> as.vector(Recall(n - 1, r - 1,
> 1:(n - 1))), frame = frame)
> if(!exists(i2))
> assign(i2, as.vector(Recall(n - 1, r,
> 1:(n - 1))), frame = frame)
> rbind(cbind(v[1],
> matrix(v[-1][get(i1)], ncol=r-1)),
> matrix(v[-1][get(i2)], ncol=r))
> } else matrix(v[1:n], ncol = 1)
> }
>
> nextSubset <- function(n, r, v = 1:n) {
> nam <- paste(n, r)
> if(exists(nam, frame = 0)) {
> ind <- get(nam)
> s <- seq(along = ind)[ind]
> if(s[1] == n - r + 1) {
> remove(nam, frame = 0)
> return(NULL)
> }
> k <- max(s[s < (n-r+1):n])
> t <- sum(ind[k:n])
> ind[k:n] <- c(F, rep(T, t),
> rep(F, n - k - t))
> } else {
> ind <- c(rep(T, r), rep(F, n-r))
> }
> assign(nam, ind, frame = 0)
> v[1:n][ind]
> }
>
>
> Pam
>
> ###################################
> Pam Goodman
> Biometrician
> Northwest Indian Fisheries Commission
> 6730 Martin Way E
> Olympia, WA 98516-5540
> Phone: (360)438-1181 ext 380
> Fax: (360)753-8659
> http://www.nwifc.wa.gov/
> ###################################
> -----------------------------------------------------------------------
> 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
>
-----------------------------------------------------------------------
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
|