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

[Prof Brian Ripley]

> From: Laurent Gautier <Laurent.Gautier at lionbioscience.com>
> Date: Tue, 4 Apr 2000 14:51:12 +0200
> 
> 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)...
> 
> # 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


prod(dlnorm(x,meanlog=a,sdlog=b)) is a single number.  Indeed, 
I suspect that dlnorm(x,meanlog=a,sdlog=b) is not what you expected:
 [1] 0.06006905 0.03601262 0.01952866 0.02211013 0.01786005 0.01354265
 [7] 0.01106827 0.00982311 0.02555533 0.01979742 0.01535234 0.01206240

even though x is of length 10, as a and b are of length 12.  You could use
something like

A <- expand.grid(x=x, meanlog=a, sdlog=b)
res <- do.call("dlnorm", A)
dim(res) <- c(length(x), length(a), length(b))
apply(res, c(2,3), prod)
             [,1]         [,2]         [,3]         [,4]
[1,] 1.141655e-12 4.952102e-14 5.144870e-15 8.780997e-16
[2,] 2.823226e-13 2.661372e-14 3.628082e-15 7.021947e-16
[3,] 9.448602e-15 5.880062e-15 1.551789e-15 4.077529e-16

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

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