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 [EMAIL PROTECTED]
>
> /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@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:[EMAIL PROTECTED] UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901
______________________________________________
R-help@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.
