I am debugging some code and found a function that returns and error most of the time. I have issolated the problem to when it raises an argument to the 1/3 power but for the life of me I can not figure out why it is not working. i have gone through the FAQs and have found nothing (not to say I might have missed something). I am highly embarrased with my inability find the problem (my face is bright red!)

I am using R 1.8.1 on win xp.

Below is the resluts of my manual testing:

> Wq <- (1+k*(Qa-(k/(6*n))))^(1/3)
> Wq
[1] NaN
> Wq <- (1+k*(Qa-(k/(6*n))))
> Wq
[1] -10.72508
> Wq <- Wq^(1/3)
> Wq
[1] NaN
> z <- -10.72508^(1/3)
> z
[1] -2.205296


as you can see if Wq = -10.72508,  Wq^(1/3) is NaN but
z<- -10.72508^(1/3) returns a number.

If someone can explain what I am doing incorrectly I would be most grateful.

Below is the complete code for the function for reference:

HallBoot <- function (x) {
#statisics from the original data
psd <- sqrt((sum((x-mean(x))^2)/length(x)))
pmean <- mean(x)
n <- length(x)
k <- ((sum((x-mean(x))^3)/(length(x)*psd^3)))
#the calulation of Q
Qstat <- function (x,mean){
	nb = length(x)
	bmean <- mean(x)
	Sb <- sqrt((sum((x-mean(x))^2)/nb))
 	Kb <- ((sum((x-mean(x))^3)/(nb*psd^3)))
	W <- (bmean-mean)/Sb
	Q <- W + (Kb * W^2)/3 + (Kb^2)*(W^3)/27 + Kb/(6*nb)
	}

#Bootsrap 10000 values for Q
Q <- bootstrap(x,10000,Qstat,pmean) 

#Determine the quantile of Q to use
Qa <- quantile(Q$thetastar, 0.05, na.rm = TRUE, names = FALSE)

#Calulate Wq*******Here is the code tha does not work***************
Wq <- (3/k)*((1+k*(Qa-(k/(6*n))))^(1/3)-1)

#Halls bootstrap 
UCL <- pmean-(Wq*psd)
UCL
}
 




Michael J. Bock, PhD.
ARCADIS
24 Preble St. Suite 100
Portland, ME 04101
207.828.0046
fax 207.828.0062




	[[alternative HTML version deleted]]


