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