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

Date: Fri 28 Oct 2005 - 02:32:11 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 05:33:58 2005

Date: Fri 28 Oct 2005 - 02:32:11 EST

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 Received on Fri Oct 28 05:33:58 2005

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