[R] case study of efficient R coding

[Jason Liao]

The following are four ways I coded for the variance of multinomial
distribution and their CPU time. It is a bit surprising that the 3rd
method beats the 4th one substantially.


 var.multinomial1 = function(n, pi)
+   {
+      p = length(pi);
+      var = array(0, c(p,p));
+      
+      for(i in 1:p)
+        for(j in 1:p)
+        {
+           if(i==j) var[i,j] = n*pi[i]*(1-pi[1])
+   else var[i,j] = - n*pi[i]*pi[j];
+        }
+    }
> 
>    var.multinomial2 = function(n, pi)
+    {
+      n*( -pi%*%t(pi) + diag(pi) )
+    }
>    
>            var.multinomial3 = function(n, pi)
+            {
+               var = -pi%*%t(pi);
+               var[row(var)==col(var)] = pi*(1-pi); 
+               n*var;
+            }
>    
>    var.multinomial4 = function(n, pi)
+            {
+               var = -pi%*%t(pi);
+               var[seq(1, by=(p+1), length=p)] = pi*(1-pi); 
+               n*var;
+            }
>    
>    n = 100;
>    pi = exp( rnorm(10));
>    pi = pi/sum(pi);
>    n.simu = 1000;
>    
>    system.time(for(i in 1:n.simu) var = var.multinomial1(n, pi))
[1]  4.50  0.04 99.73    NA    NA
>     system.time(for(i in 1:n.simu) var = var.multinomial2(n, pi))
[1]  0.33  0.00 10.74    NA    NA
>      system.time(for(i in 1:n.simu) var = var.multinomial3(n, pi))
[1] 0.14 0.00 4.17   NA   NA
>       system.time(for(i in 1:n.simu) var = var.multinomial4(n, pi))
[1] 0.21 0.00 6.47   NA   NA
>    
>      
>    
But first, does anyone know a vectorized way of coding (summing the
outter product of rows)

p = 10;
x = rnorm(p*p);
dim(x) = c(p.p);

var = array(0, c(p,p));
for(i in 1:p) var = var + x[i,]%*% t(x[i,]);

=====
Jason G. Liao, Ph.D.
Division of Biometrics
University of Medicine and Dentistry of New Jersey
335 George Street, Suite 2200
New Brunswick, NJ 08903-2688
phone (732) 235-8611, fax (732) 235-9777
http://www.geocities.com/jg_liao

__________________________________________________

Faith Hill - Exclusive Performances, Videos & More

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._