[R] outer-question

[Rau, Roland]

Dear all,

a big thanks to Thomas Lumley, James Holtman and Tony Plate for their
answers. They all pointed in the same direction => I need a vectorized
function to be applied. Hence, I will try to work with a 'wrapper'
function as described in the FAQ.

Thanks again,
Roland


> -----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}}
> >
> > ______________________________________________
> > 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
> >
> 
> Thomas Lumley			Assoc. Professor, Biostatistics
> tlumley at u.washington.edu	University of Washington, Seattle
> 

+++++
This mail has been sent through the MPI for Demographic Rese...{{dropped}}