From: Thomas Lumley <tlumley_at_u.washington.edu>

Date: Fri 28 Oct 2005 - 07:39:21 EST

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Fri Oct 28 08:02:38 2005

Date: Fri 28 Oct 2005 - 07:39:21 EST

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@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@u.washington.edu University of Washington, Seattle ______________________________________________R-help@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 Received on Fri Oct 28 08:02:38 2005

*
This archive was generated by hypermail 2.1.8
: Fri 28 Oct 2005 - 09:48:26 EST
*