thr3ads.net - R help - [R] outer-question [Oct 2005]

If this information is useful, please help other people find it:
Share via:

Rau, Roland

2005-Oct-27 16:32 UTC

[R] outer-question

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}}

Tony Plate

2005-Oct-27 17:19 UTC

head link

[R] outer-question

It looks like you didn't vectorize the function you gave "outer"
in your
longer example.

Consider your short example with a diagnostic printout:

 > a <- 1:3
 > b <- 1:4
 > f <- function(a,b,d) {
+     cat("In f:", length(a), length(b), "\n")
+     return(a*b+(sum(d)))
+ }
 > additional <- runif(100)
 > outer(X=a, Y=b, FUN=f, d=additional)
In f: 12 12
          [,1]     [,2]     [,3]     [,4]
[1,] 53.61985 54.61985 55.61985 56.61985
[2,] 54.61985 56.61985 58.61985 60.61985
[3,] 55.61985 58.61985 61.61985 64.61985
 >

Note that "f" is called only once, with vectors for "a" and
"b".

-- Tony Plate

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

2005-Oct-27 21:39 UTC

head link

[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

Rau, Roland

2005-Oct-28 07:44 UTC

head link

[R] outer-question

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}}

Reasonably Related Threads

Search for more reasonably related threads

R help - Oct 2005 - outer-question

[R] outer-question

[R] outer-question

[R] outer-question

[R] outer-question

Reasonably Related Threads