0

我在 OSX 中执行我的 R 代码时遇到问题。那是我的代码：

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)

当我初始化 x,y,q,p 并使用日志时，它可以工作。

但是为什么会出现这些错误，但是为什么？

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
4

3 回答 3

4

这是另一种方法。如果您减少需要迭代地做某事的次数，您的代码通常应该更快。具体来说，您的所有runif(1,0,1)调用都可以替换为一个大的runif()值向量，然后根据该向量对向量进行子集化。

我以@Mark Miller 的函数为起点，进行了以下修改。请注意，如果过采样器保留前一组随机数中的良好值并且只填充直到达到，这可以进一步改进n，但无论如何这都非常快。对于速度比较，我逐字记录了他的代码并将其包装在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")

和结果

           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

从新的 R 会话开始，复制和粘贴上面的代码不会返回错误。这是一个例子：

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

在此处输入图像描述

在此处输入图像描述

于 2012-04-14T12:29:41.893 回答
2

这是做你想做的吗？一世：

  1. 添加向量来存储 x、y、q 和 p。
  2. 将 ln 更改为 log
  3. 添加了一个 plot(a,b) 语句。
  4. 出于调试目的将 20000 更改为 20。

我没有 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)
于 2012-04-14T08:42:48.383 回答
2

您不能从 Windows 上一个全新的空工作区开始。' x' 对象必须已经存在，否则你也会在那里得到错误。在 Windows 和 OSX 中执行ls()并查看是否有 ' x'。我会把钱放在 Windows 上而不是 OSX 上。

于 2012-04-14T09:46:18.893 回答