[R] Question about the proper use of outer()

[Peter Dalgaard BSA]

Laurent Gautier <Laurent.Gautier at lionbioscience.com> writes:

> DeaR all,
> 
> 
> I do not have a clue with is the following NOT working like I expect to
> do...
> (and I cannot find any answer at CRAN)...

We actually had this recently... Exchange between Doug Bates and Brian
Ripley IIRC. 

> # This one is for my sample
> > x _ seq(3,10)
> # This two for parameters
> > a _ seq(2,4)
> > b _ seq(2,5)
> #  This one for the likelihood of a sample
> >f _ function(a,b) {
> +  prod(dlnorm(x,meanlog=a,sdlog=b))
> + }
> >  outer(a,b,f)
>              [,1]         [,2]         [,3]         [,4]
> [1,] 2.301348e-21 2.301348e-21 2.301348e-21 2.301348e-21
> [2,] 2.301348e-21 2.301348e-21 2.301348e-21 2.301348e-21
> [3,] 2.301348e-21 2.301348e-21 2.301348e-21 2.301348e-21
> 
> # interestingly, the following seems to work:
> 
> > f(a[1],b[1])
> [1] 1.141655e-12
> > f(a[1],b[2])
> [1] 4.952102e-14
> 

The problem is that your f doesn't vectorize in a and b. It gets
called with 

a = c(2,3,4,2,3,4,2,3,4,2,3,4)
b = c(2,2,2,3,3,3,4,4,4,5,5,5)

(or vice versa, you get the picture) Then it gets to the
dlnorm(x,meanlog=a,sdlog=b) and recycles x to get length 12 for all
vectors and takes the product of everything. The easy, although not
the most efficient, way out is to vectorize explicitly:

g<-function(a,b) sapply(seq(along=a), function(i)f(a[i],b[i]))
outer(a,b,g)


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