[R] vectorize a function

Christos Hatzis christos at nuverabio.com
Fri Jun 22 16:17:24 CEST 2007


How about:

sum(sapply(unique(a), function(x) {b <- which(a==x); sum(M[b, b])}))

HTH
-Christos

Christos Hatzis, Ph.D.
Nuvera Biosciences, Inc.
400 West Cummings Park
Suite 5350
Woburn, MA 01801
Tel: 781-938-3830
www.nuverabio.com
 

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Robin Hankin
> Sent: Friday, June 22, 2007 9:28 AM
> To: RHelp help
> Subject: [R] vectorize a function
> 
> Hello everyone
> 
> suppose I have an integer vector "a" of length "n" and a 
> symmetric matrix "M" of size n-by-n.
> 
> Vector "a" describes a partition of a set of "n" elements and 
> matrix M describes a penalty function: row i column j 
> represents the penalty if element i and element j are in the 
> same partition.
> 
> Toy example follows; the real case is much larger and I need 
> to evaluate my penalty function many times.
> 
> If a <- c(1,1,2,1,3)  then elements 1,2,4 are in the same 
> partition; element 3 is in a partition on its own and element 
> 5 is in a partition on its own.
> 
> The total penalty  can be described by the following (ugly)
> function:
> 
> f <- function(a,M){
>    out <- 0
>    for(i in unique(a)){
>      out <- out + sum(M[which(a==i),which(a==i)])
>    }
>    return(out)
> }
> 
> 
> so with
> 
> M <- matrix(rpois(25,3),5,5)
> M <- M+t(M)
> diag(M) <- 0
> a <- c(1,2,1,1,3)
> 
> f(a,M) gives the total penalty.
> 
> 
> QUESTION:  how to rewrite f() so that it has no loop?
> 
> 
> 
> 
> 
> 
> --
> Robin Hankin
> Uncertainty Analyst
> National Oceanography Centre, Southampton European Way, 
> Southampton SO14 3ZH, UK
>   tel  023-8059-7743
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
>



More information about the R-help mailing list