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:[EMAIL PROTECTED] > Sent: Thursday, October 27, 2005 11:39 PM > To: Rau, Roland > Cc: [email protected] > 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}} > > > > ______________________________________________ > > [email protected] 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 > [EMAIL PROTECTED] University of Washington, Seattle > +++++ This mail has been sent through the MPI for Demographic Rese...{{dropped}} ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
