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