蒙特卡罗模拟代码：在 R 中生成给定大小的样本答案

【问题标题】：Code for Monte Carlo simulation: generate samples of given size in R蒙特卡罗模拟代码：在 R 中生成给定大小的样本
【发布时间】：2013-04-25 17:46:18
【问题描述】：

我首先使用以下代码生成 500 个在 0 和 1 之间均匀分布的随机数的样本：

set.seed(1234)
X<-runif(500, min=0, max=1)

现在，我需要编写一个伪代码，为 MC 模拟生成 10000 个 N=500 的样本，计算我新创建的 X 的平均值，并将迭代次数和平均值存储在结果对象中。我从来没有尝试过，到目前为止我有这个：

n.iter <-(10000*500)
results <- matrix (0, n.iter, 4)

最后，一旦完成，我将运行它，然后获取累积样本均值的中值、均值和最小/最大值，并将它们保存到名为 MC.table 的数据帧中。（还要注意，上面，我不知道为什么矩阵代码中有一个“4”——我正在研究前面的例子）。任何建议或帮助将不胜感激。

编辑：我有一个可能有用的例子，但我不太明白它是怎么回事，所以请详细说明它的有效性：

Ni <- 10000
n <- 500
c <- 0

for (i in n){
for (j in 1:Ni){
c <- c+ 1
d <- data.frame (x= , y= )
results [c,1] <- c
results [c,2] <- j
results [c,3] <- i
results [c,4] <- something( d$x, d$y)
rm (d) } }

如果你能花时间解释这意味着什么，那对我有很大帮助！谢谢！

【问题讨论】：

您需要很好地掌握 R 中的数据结构和“for 循环”。我从 Matloff 的 R for Programmers 中学到了很多东西：heather.cs.ucdavis.edu/~matloff/R/RProg.pdf
谢谢，我将不得不对此进行更多研究。我是一个认真的 R 初学者，我的教授在对材料进行了一次讲座并使用 SAS 参考材料后让我们进行 MC 模拟。如果我尝试过，我不会再迷失了！
编辑您的问题并添加详细信息。
我认为你在正确的轨道上。通往同一个目的地的道路有很多。如果与下面的代码进行比较（出于学习 R 的目的），我将内部 for 循环替换为函数“sample”。

标签： r montecarlo

【解决方案1】：

您可以尝试使用data.table，这是一个可以使用install.packages("data.table") 安装的软件包。安装后，您将运行类似...

> require(data.table)
> dt <- data.table(x=runif(500*10000),iter=rep(1:500,each=10000))
                  # x iter
      # 1: 0.48293196    1
      # 2: 0.61935416    1
      # 3: 0.99831614    1
      # 4: 0.26944687    1
      # 5: 0.38027524    1
     # ---                
# 4999996: 0.11314160  500
# 4999997: 0.07958396  500
# 4999998: 0.97690312  500
# 4999999: 0.81670765  500
# 5000000: 0.62934609  500
> summaries <- dt[,list(mean=mean(x),median=median(x)),by=iter]
     # iter      mean    median
  # 1:    1 0.5005310 0.5026592
  # 2:    2 0.4971551 0.4950034
  # 3:    3 0.4977677 0.4985360
  # 4:    4 0.5034727 0.5052344
  # 5:    5 0.4999848 0.4971214
 # ---                         
# 496:  496 0.5013314 0.5048186
# 497:  497 0.4955447 0.4941715
# 498:  498 0.4983971 0.4910115
# 499:  499 0.5000382 0.4997024
# 500:  500 0.5009614 0.4988237
> min_o_means <- min(summaries$mean)
# [1] 0.4914826

我认为语法相当简单。您可能想使用?（例如?rep）查找一些函数。以 # 开头的行只是显示生成的对象。在 data.tables 中，: 左边的数字只是行号，--- 表示在显示中被跳过的行。

【讨论】：

很好，但我认为这对于初学者来说有点过头了。 :) 你也可以添加for 和apply 解决方案。
说实话，我不知道标准的应用解决方案，也不擅长编写优雅的 for 循环。 :) 我希望其他人会为 OP 写这些。在我开始使用 data.table 之前，我的代码一团糟。
data.tables 的好方法我也遇到了应用问题。
@RomanLuštrik 所以我刚刚了解了社区维基。将我的答案更改为一个好主意，以便有人可以添加您建议的部分吗？

【解决方案2】：

我想我会给出的答案真的取决于你是想学习伪代码还是你想学习“R”式的方法。这个答案是我向想要学习如何使用 R 的人推荐的。

首先，我将创建一个包含 N 列和 10000 行的矩阵。当我们提前腾出空间让数字进入时，R 会很感激。

X=matrix(NA,nrow=10000,ncol=500)

您知道如何为一行生成 500 个随机变量。

runif(500,0,1)

现在你需要弄清楚如何让它发生 10000 次并将每个分配给 X。也许 for 循环会起作用。

for(i in 1:10000) X[i,]=runif(500,0,1)

然后你需要弄清楚如何获取每一行的摘要。一个可能有帮助的功能是rowMeans()。查看它的帮助页面，然后尝试获取表格 X 每一行的平均值

获取每次迭代的均值

rowMeans(X)

然后为了了解这些数字是什么样的，我可能倾向于

plot(rowMeans(X))

【讨论】：

【解决方案3】：

我认为您正在描述一个简单的引导程序。最终，您可能想要使用启动功能。但在你理解机制之前，我觉得循环是要走的路。这应该让你开始：

test<-function(
    seed=1234,
    sample.size=500,
    resample.number=1000,
    alpha=0.05
    )
    {

        #initialize original sample
        original.sample<-runif(sample.size, min=0, max=1)   

        #initialize data.frame
        resample.results<-data.frame("Run.Number"=NULL,"mean"=NULL)
        for(counter in 1:resample.number){
            temp<-sample(original.sample, size=length(original.sample), replace = TRUE)
            temp.mean<-mean(temp)
            temp.table.row<-data.frame("Run.Number"=counter,"mean"=temp.mean)
            resample.results<-rbind(resample.results,temp.table.row)
        }
        resample.results<-resample.results[with(resample.results, order(mean)), ]

        #for the mean information
        lowerCI.row<-resample.number*alpha/2
        upplerCI.row<-resample.number*(1-(alpha/2))
        median.row<-resample.number/2

        #for the mean information
        median<-resample.results$mean[median.row]
        lowerCI<-resample.results$mean[lowerCI.row]
        upperCI<-resample.results$mean[upplerCI.row]

        #for the position information
        median.run<-resample.results$Run.Number[median.row]
        lowerCI.run<-resample.results$Run.Number[lowerCI.row]
        upperCI.run<-resample.results$Run.Number[upplerCI.row]

        mc.table<-data.frame("median"=NULL,"lowerCI"=NULL,"upperCI"=NULL)
        values<-data.frame(median,lowerCI,upperCI)
        #as.numeric because R doesn't like to mix data types
        runs<-as.numeric(data.frame(median.run,lowerCI.run,upperCI.run))
        mc.table<-rbind(mc.table,values)
        mc.table<-rbind(mc.table,runs)

        print(mc.table)
    }

对数据重新采样后，您会找到均值。然后您订购所有重新采样的方法。该列表的中间是中位数。并且，例如，对于 10000 个重采样，第 250 个有序重采样均值将是您较低的 95% CI。虽然我在这里没有这样做，但最小值将在位置 1，最大值将在位置 10000。降低重采样数时要小心：我计算位置的方式可能会变成十进制值，这会混淆R.

顺便说一句，我把它放在函数形式中。如果您喜欢逐行处理，只需确保运行 function() 和以下主要 {}

之间的所有行

【讨论】：