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]]