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]]
Well, your example can be simplified to> -1^(1/3)[1] -1> (-1)^(1/3)[1] NaN Do you see the light? [think "precedence rules"] Regards, Martin
"Bock, Michael" <MBock at arcadis-us.com> writes:> > 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.Hint:> -2^(1/2)[1] -1.414214> sqrt(-2)[1] NaN Warning message: NaNs produced in: sqrt(-2)> sqrt(-2+0i)[1] 0+1.414214i Powers of negative numbers are not defined for non-integral exponents. You need something like cuberoot <- function(x)sign(x)*abs(x)^(1/3) -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
I see the light! I have gotten response re precedence : i.e (-2)^2 = 4 and -2^2 = -4 (^ evaluated before -) I understood precedence and still could not see the problem but as Peter states: Powers of negative numbers are not defined for non-integral exponents. You need something like cuberoot <- function(x)sign(x)*abs(x)^(1/3) I though you could take the cubed root of negative numbers (basic algebra) but perhaps this does not hold true for computer math due to rounding of the exponents. I am on my way! . Thanks 10^6 Mike -----Original Message----- From: Peter Dalgaard [mailto:p.dalgaard at biostat.ku.dk] Sent: Tuesday, April 27, 2004 12:17 PM To: Bock, Michael Cc: R-help at stat.math.ethz.ch Subject: Re: [R] Problems raised to 1/3 power and NaN "Bock, Michael" <MBock at arcadis-us.com> writes:> > 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.Hint:> -2^(1/2)[1] -1.414214> sqrt(-2)[1] NaN Warning message: NaNs produced in: sqrt(-2)> sqrt(-2+0i)[1] 0+1.414214i Powers of negative numbers are not defined for non-integral exponents. You need something like cuberoot <- function(x)sign(x)*abs(x)^(1/3) -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
On Tue, 27 Apr 2004, Bock, Michael wrote:> > I though you could take the cubed root of negative numbers (basic > algebra) but perhaps this does not hold true for computer math due to > rounding of the exponents. >Yes. All the representable floating point numbers are fractions whose denominator is a power of two, and so none of them give exactly real-valued roots of negative numbers. You could use complex numbers, but then you have the problem of getting the correct cube root.> (-1+0i)^(1/3)[1] 0.5+0.8660254i -thomas