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

[Sergey Goriatchev]

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?

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
       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)
}
#======================================================