outer-question

-----Original Message-----
From: Thomas Lumley [mailto:tlumley at u.washington.edu] 
Sent: Thursday, October 27, 2005 11:39 PM
To: Rau, Roland
Cc: r-help at stat.math.ethz.ch
Subject: Re: [R] outer-question


You want FAQ 7.17 Why does outer() behave strangely with my function?

 	-thomas

On Thu, 27 Oct 2005, 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}}