From: Rau, Roland <Rau_at_demogr.mpg.de>

Date: Fri 28 Oct 2005 - 17:44:25 EST

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 17:50:01 2005

Date: Fri 28 Oct 2005 - 17:44:25 EST

> -----Original Message-----

*> From: Thomas Lumley [mailto:tlumley@u.washington.edu]
**> Sent: Thursday, October 27, 2005 11:39 PM
**> To: Rau, Roland
**> Cc: r-help@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@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
**>
*

+++++

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 Received on Fri Oct 28 17:50:01 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:40:51 EST
*