[Rd] [Fwd: Re: [R] a strange problem with integrate()]

[Patrick Burns]

When I saw the subject of the original message on
R-help, I was 95% confident that I knew the answer
(before I had seen the question).

This made me think that perhaps for some functions
there should be a 'Troubleshooting' section in the help
file.

The current help file for 'integrate' does say, as Sundar
points out, what the requirements are.  However, I
think more people would solve the problem more quickly
on their own if there were a troubleshooting section.

Most functions aren't abused in predictable ways, but
a few are.  Another that springs immediately to mind
is 'read.table'.

Patrick Burns
patrick at burns-stat.com
+44 (0)20 8525 0696
http://www.burns-stat.com
(home of S Poetry and "A Guide for the Unwilling S User")

-------- Original Message --------
Subject: 	Re: [R] a strange problem with integrate()
Date: 	Wed, 01 Mar 2006 11:44:33 -0600
From: 	Sundar Dorai-Raj <sundar.dorai-raj at pdf.com>
Organization: 	PDF Solutions, Inc.
To: 	vito muggeo <vmuggeo at dssm.unipa.it>
CC: 	r-help at stat.math.ethz.ch <r-help at stat.math.ethz.ch>
References: 	<4405D98F.6030709 at dssm.unipa.it>



vito muggeo wrote:
> Dear all,
> I am stuck on the following problem with integrate(). I have been out of 
> luck using RSiteSearch()..
> 
> My function is
> 
> g2<-function(b,theta,xi,yi,sigma2){
>        xi<-cbind(1,xi)
>        eta<-drop(xi%*%theta)
>        num<-exp((eta + rep(b,length(eta)))*yi)
>        den<- 1 + exp(eta + rep(b,length(eta)))
>        result=(num/den)*exp((-b^2)/sigma2)/sqrt(2*pi*sigma2)
>        r=prod(result)
>        return(r)
>        }
> 
> And I am interested in evaluating the simple integral, but:
> 
>  > integrate(g2,-2,2,theta=c(-2,3),xi=c(1,2,5,6),yi=c(1,0,1,1),sigma2=1)
> Error in integrate(g2, -2, 2, theta = c(-2, 3), xi = c(1, 2, 5, 6), yi = 
> c(1,  :
>          evaluation of function gave a result of wrong length
>  >
> 
> I have checked the integrand function
> 
>  > valori<-seq(-2,2,l=30)
>  > risvalori<-NULL
>  > for(k in valori) 
> risvalori[length(risvalori)+1]<-g2(k,theta=c(-2,3),xi=c(1,2,5,6),yi=c(1,0,1,1),sigma2=1)
>  > plot(valori, risvalori)
> 
> And the integral exists..
> 
> Please, any comment is coming?
> 
> many thanks,
> vito
> 

Please (re-)read ?integrate:

  f: an R function taking a numeric first argument and returning a
           numeric vector of the same length.  Returning a non-finite
           element will generate an error.

Note the "returning a numeric vector of the *same* length." Your 
function returns "prod(r)" which is not the same length as "b".

Some style issues (and I state these as diplomatically as is possible in 
e-mail):

a. Don't mix "<-" with "=" for assignment in your scripts.
b. Use more spaces and consistent indenting.

Here's what my code looks like:

g2 <- function(b, theta, xi, yi, sigma2) {
   xi <- cbind(1, xi)
   eta <- drop(xi %*% theta)
   num <- exp((eta + rep(b, length(eta))) * yi)
   den <- 1 + exp(eta + rep(b, length(eta)))
   result <- (num/den) * exp((-b2)/sigma2)/sqrt(2 * pi * sigma2)
   r <- prod(result)
   r
}

After reformatting your code I saw your problem immediately without 
having executing a single line.

HTH,

--sundar

______________________________________________
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