CODEWITHSUNDEEP
macosrstatisticsplot
jeffrey
jeffrey

Reputation: 2096

Code works fine on win7, but there are problems on OSX

i have problems to perform my R-Code in OSX. That's my code:

i <- 1
while (i <= 20000) {
  repeat{
    z1=((runif(1,0,1)*2)-1)
    z2=((runif(1,0,1)*2)-1)
    h=z1**2+z2**2
    if((h > 0) && (h <= 1)){break}
   }
  x[i] <- z1
  y[i] <- z2
  q[i] <- h

  i <- i + 1
} 

j <- 1
while (j <= 20000) {
  h=sqrt((-2*ln(q[j]))/q[j])
  p[j] <- h
  j <- j + 1
}

a=x*p
b=y*p
points(a,b, pch=c(20,20),col=c("dark green","red"),cex=0.6)

When I initialize x,y,q,p and use the log, it works.

But why are there those errors, but why?

error in x[i] <- z1: object 'x' not fund
error: no function for "ln" fund
error: object 'x' not fund
error: object 'y' not fund

Upvotes: 0

Views: 172

Answers (3)

Chase
Chase

Reputation: 69201

Here's another approach. Your code should generally be faster if you reduce the number of times you need to do things iteratively. Specifically, all of your runif(1,0,1) calls can be replaced by one big vector of runif() values, then subset the vector based on that.

I used @Mark Miller's function as the starting point and made the following modifications. Note, this can be improved further if the oversampler kept the good values from the previous set of random numbers and only filled in until n was reached, but this is pretty fast regardless. For the speed comparisons, I took his code verbatim and wrapped it in fun2 <- function() {...}

fun1 <- function(n, oversample = 1.50){
  #oversample
  over <- ceiling(n * oversample)
  goodvars <- NA
  while (length(goodvars) < n){
    z1 <- runif(over,-1,1)
    z2 <- runif(over,-1,1)
    h <- z1^2 + z2^2  
    goodvars <- which(h > 0 & h < 1)
  }
  goodvars <- goodvars[1:n]
  x <- z1[goodvars]
  y <- z2[goodvars]
  q <- h[goodvars]
  p <- sqrt((-2 * log(q)) / q)
  a <- x * p
  b <- y * p
  return(cbind(a,b))
 }

##Mark's code put into a function
fun2 <- function() {
  i <- 1

  x <- rep(NA, 20)
  y <- rep(NA, 20)
  q <- rep(NA, 20)
  p <- rep(NA, 20)

  while (i <= 20) {

    repeat{
      z1=((runif(1,0,1)*2)-1)
      z2=((runif(1,0,1)*2)-1)
      h=z1**2+z2**2
      if((h > 0) & (h <= 1)){break}
    }
    x[i] <- z1
    y[i] <- z2
    q[i] <- h

    i <- i + 1
  } 

  j <- 1
  while (j <= 20) {

    h=sqrt((-2*log(q[j]))/q[j])

    p[j] <- h

    j <- j + 1
  }

  a=x*p
  b=y*p
}

#Do some speed checking with rbenchmark. Also checkout compiler package for some free speed
library(compiler)
library(rbenchmark)
#Compile functions to see improvements
cfun1 <- cmpfun(fun1)
cfun2 <- cmpfun(fun2)
#run benchmark tests
benchmark(fun1(n = 20), fun2(), cfun1(n = 20), cfun2(),
          replications = 1000,
          columns=c("test", "elapsed", "relative"),
          order = "elapsed")

And the results

           test elapsed  relative
3 cfun1(n = 20)   0.042  1.000000
1  fun1(n = 20)   0.055  1.309524
4       cfun2()   0.407  9.690476
2        fun2()   0.882 21.000000

Starting with a new R session, copying and pasting the code above does not return an error. Here's an example:

test <- fun1(n = 1000)
plot(test)

enter image description here

enter image description here

Upvotes: 4

Spacedman
Spacedman

Reputation: 94202

You can't be starting from a fresh empty workspace on Windows. The 'x' object must already exist, or you'd get the error there too. Do ls() in Windows and OSX and see if there is an 'x'. I'd put money on there being one in Windows but not OSX.

Upvotes: 2

Mark Miller
Mark Miller

Reputation: 13113

Does this do what you want? I:

  1. added vectors to store x, y, q, and p.
  2. changed ln to log
  3. added a plot(a,b) statement.
  4. changed 20000 to 20 for debugging purposes.

I do not have a Mac.

i <- 1

x <- rep(NA, 20)
y <- rep(NA, 20)
q <- rep(NA, 20)
p <- rep(NA, 20)

while (i <= 20) {

repeat{
    z1=((runif(1,0,1)*2)-1)
    z2=((runif(1,0,1)*2)-1)
    h=z1**2+z2**2
if((h > 0) & (h <= 1)){break}
 }
x[i] <- z1
y[i] <- z2
q[i] <- h

i <- i + 1
} 

j <- 1
while (j <= 20) {

h=sqrt((-2*log(q[j]))/q[j])

p[j] <- h

j <- j + 1
}

a=x*p
b=y*p
plot(a,b)
points(a,b, pch=c(20,20),col=c("dark green","red"),cex=0.6)

Upvotes: 2

Related Questions