[R] outer-question

[Tony Plate]

It looks like you didn't vectorize the function you gave "outer" in your 
longer example.

Consider your short example with a diagnostic printout:

 > a <- 1:3
 > b <- 1:4
 > f <- function(a,b,d) {
+     cat("In f:", length(a), length(b), "\n")
+     return(a*b+(sum(d)))
+ }
 > additional <- runif(100)
 > outer(X=a, Y=b, FUN=f, d=additional)
In f: 12 12
          [,1]     [,2]     [,3]     [,4]
[1,] 53.61985 54.61985 55.61985 56.61985
[2,] 54.61985 56.61985 58.61985 60.61985
[3,] 55.61985 58.61985 61.61985 64.61985
 >

Note that "f" is called only once, with vectors for "a" and "b".

-- Tony Plate

Rau, Roland wrote:
> Dear all,
> 
> This is a rather lengthy message, but I don't know what I made wrong in
> my real example since the simple code works.
> I have two variables a, b and a function f for which I would like to
> calculate all possible combinations of the values of a and b.
> If f is multiplication, I would simply do:
> 
> a <- 1:5
> b <- 1:5
> outer(a,b)
> 
> ## A bit more complicated is this:
> f <- function(a,b,d) {
> 	return(a*b+(sum(d)))
> }
> additional <- runif(100)
> outer(X=a, Y=b, FUN=f, d=additional)
> 
> ## So far so good. But now my real example. I would like to plot the
> ## log-likelihood surface for two parameters alpha and beta of 
> ## a Gompertz distribution with given data
> 
> ### I have a function to generate random-numbers from a
> Gompertz-Distribution
> ### (using the 'inversion method')
> 
> random.gomp <- function(n, alpha, beta) {
>             return( (log(1-(beta/alpha*log(1-runif(n)))))/beta)
> }
> 
> ## Now I generate some 'lifetimes'
> no.people <- 1000
> al <- 0.1
> bet <- 0.1
> lifetimes <- random.gomp(n=no.people, alpha=al, beta=bet)
> 
> ### Since I neither have censoring nor truncation in this simple case,
> ### the log-likelihood should be simply the sum of the log of the
> ### the densities (following the parametrization of Klein/Moeschberger
> ### Survival Analysis, p. 38)
> 
> loggomp <- function(alphas, betas, timep) {
>   return(sum(log(alphas) + betas*timep + (alphas/betas *
> (1-exp(betas*timep)))))
> }
> 
> ### Now I thought I could obtain a matrix of the log-likelihood surface
> ### by specifying possible values for alpha and beta with the given
> data.
> ### I was able to produce this matrix with two for-loops. But I thought
> ### I could use also 'outer' in this case.
> ### This is what I tried:
> 
> possible.alphas <- seq(from=0.05, to=0.15, length=30)
> possible.betas <- seq(from=0.05, to=0.15, length=30)
> 
> outer(X=possible.alphas, Y=possible.betas, FUN=loggomp, timep=lifetimes)
> 
> ### But the result is:
> 
>>outer(X=possible.alphas, Y=possible.betas, FUN=loggomp,
> 
> timep=lifetimes)
> Error in outer(X = possible.alphas, Y = possible.betas, FUN = loggomp,
> : 
>         dim<- : dims [product 900] do not match the length of object [1]
> In addition: Warning messages:
> ...
> 
> ### Can somebody give me some hint where the problem is?
> ### I checked my definition of 'loggomp' but I thought this looks fine:
> loggomp(alphas=possible.alphas[1], betas=possible.betas[1],
> timep=lifetimes)
> loggomp(alphas=possible.alphas[4], betas=possible.betas[10],
> timep=lifetimes)
> loggomp(alphas=possible.alphas[3], betas=possible.betas[11],
> timep=lifetimes)               
> 
> 
> ### I'd appreciate any kind of advice.               
> ### Thanks a lot in advance.
> ### Roland
>                
> 
> +++++
> This mail has been sent through the MPI for Demographic Rese...{{dropped}}
> 
> ______________________________________________
> 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
>