I have the following code for generating random variables in R
theta<-1
n<-1000
u<-runif(n)
H<-(-theta*log(1-u)^(1/3))
My problem is that the output I'm getting for H is NaN 1000 times as the output when I know that the answer is real. Is there something wrong with the way I defined H?
In order:
u = runif(n) # between 0 and 1
1 - u # between 0 and 1
log(1 - u) # negative real
log(1 - u) ^ (1/3) # NaN because...
From help("^") "Users are sometimes surprised by the value returned, for example why (-8)^(1/3) is NaN. For double inputs, R makes use of IEC 60559 arithmetic on all platforms, together with the C system function pow for the ^ operator. The relevant standards define the result in many corner cases. In particular, the result in the example above is mandated by the C99 standard."
So whenever you raise x to a non-integer power, the answer will be NaN if x < 0. One option, to get a negative real result, is to raise the abs(x) to the non-integer power and then make the result negative. (More generally, we multiply by sign(x) so it still works for positive x.)
sign(log(1 - u)) * (abs(log(1 - u)) ^ (1/3)) # Negative real, based on your inputs
-theta * sign(log(1 - u)) * (abs(log(1 - u)) ^ (1/3)) # Maybe what you want?