[R] Integrate() error message, I am at a loss

[Charles C. Berry]

On Wed, 12 Sep 2007, Sergey Goriatchev wrote:

> Hello!
>
> I have a problem with integrate() in my function nctspa(). Integrate
> produces an error message "evaluation of function gave a result of
> wrong length". I don't know what that means. Could anyone suggest me
> what is wrong with my function?

Sure, but you can do this yourself.

For one thing, setting

 	options( error = recover )

before you run your function will help you to see what is happening.

(viz. after the error message select 2, then figure out what f() and ff() 
are, then try ff(1))

>
> These are the examples of function calls that work OK:
> nctspa(a=1:10,n=5)
> nctspa(a=1:10, n=5, mu=2, theta=3, renorm=0)
>
> This does not work:
> nctspa(a=1:10, n=5, mu=2, theta=3, renorm=1)
>
> Many thanks in advance for your help!
> please, send reply also to sergeyg at gmail.com
>
> /Sergey
>
> Here is the function:
>
> #Computes the s.p.a. to the doubly noncentral t at value x.
> #degrees of freedom n, noncentrality parameters mu and theta.
> #==========================================================#
> nctspa <- function(a,n,mu=0,theta=0,renorm=0,rec=0){
> #Pass renorm=1 to renormalize the SPA to the pdf.
> #There is a last argument called rec. DO NOT PASS it!
>
> alpha <- mu/sqrt((1+theta/n))
> normconst <- 1
> if(renorm==1 & rec==0){
>       term1 <- integrate(nctspa, -Inf, alpha, n=n, mu=mu, theta=theta)$value
>       term2 <- integrate(nctspa, alpha, Inf, n=n, mu=mu, theta=theta)$value


Whoa! Let's read the help page for integrate:

Arguments:

        f: an R function taking a numeric first argument and returning a
           numeric vector of the same length. 
..........^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^.............


which is not what nctspa() does. It returns a list of length 2 not matter 
what the first arguement.

Did you mean something like

        	integrate(function(x,...) nctspa(x,...)$PDF, <etc> ??


HTH,

Chuck


>       normconst <- 1/(term1+term2)
> }
> cdf <- numeric()
> pdf <- cdf
> c3 <- n^2+2*n*a^2+a^4
> c2 <- (-2*mu*(a^3+n*a))/c3
> c1 <- (-n^2-n*a^2-n*theta+a^2*mu^2)/c3
> c0 <- (n*a*mu)/c3
> q <- c1/3-(c2^2)/9
> r <- 1/6*(c1*c2-3*c0)-1/27*c2^3
> b0 <- sqrt(-4*q)*cos(acos(r/sqrt(-q^3))/3)-c2/3
> t1 <- -mu+a*b0
> t2 <- -a*t1/b0/n/2
> nu <- 1/(1-2*t2)
> w <- sqrt(-mu*t1-n*log(nu)-2*theta*nu*t2)*sign(a-alpha)
> u <- sqrt((a^2+2*n*t2)*(2*n*nu^2+4*theta*nu^3)+4*n^2*b0^2)/(2*n*b0^2)
> pdf <- normconst*dnorm(w)/u
>
> nz <- (abs(t1*b0)>=1e-10)
> iz <- (abs(t1*b0)<=1e-10)
> if(any(nz)){
>       d <- numeric()
>       d[nz] <- 1/(t1[nz]*b0[nz])
>       cdf[nz] <- pnorm(w[nz])+dnorm(w[nz])*(1/w[nz]-d[nz]/u[nz])
> }
> if(any(iz)){
>       n <- sum(iz==1)
>       rez <- nctspa(c(a[iz]-1e-4, a[iz]+1e-4),n,mu,theta,0,rec+1)
>       if(rec>5){
>       cdf[iz] <- 0.5
>       warning('Too many recursions')
>       } else {
>       cdf[iz] <- 0.5*(rez$CDF[1:n]+rez$CDF[(n+1):length(rez$CDF)])
>       }
> }
> list(PDF=pdf, CDF=cdf)
> }
> #======================================================
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

Charles C. Berry                            (858) 534-2098
                                             Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901