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