Montecarlo模拟R中Stata命令的等价性

Question

我正在寻找一个非常方便的Stata命令simulate的R中的等价函数。该命令基本上允许您声明一个program (下面的示例中是reg_simulation)，然后从simulate调用这样的程序并存储所需的输出。

下面是simulate程序使用情况的Stata演示，以及我使用R复制它的尝试。

最后，我的主要问题是：是R用户运行Montecarlo模拟的方式吗？还是在结构或速度瓶颈方面遗漏了什么？提前谢谢你。

Stata实例

1. 定义reg_simulation程序。

clear all
*Define "reg_simulation" to be used later on by "simulate" command 
program reg_simulation, rclass
    *Declaring Stata version
    version 13
    *Droping all variables on memory
    drop _all
    *Set sample size (n=100)
    set obs 100
    *Simulate model
    gen x1 = rnormal()
    gen x2 = rnormal()
    gen y = 1 + 0.5 * x1 + 1.5 *x2 + rnormal()
    *Estimate OLS
    reg y x1 x2 
    *Store coefficients
    matrix B = e(b)
    return matrix betas = B 
end

1. 从reg_simulation命令调用simulate：

*Seet seed
set seed 1234
*Run the actual simulation 10 times using "reg_simulation"
simulate , reps(10) nodots: reg_simulation

1. 获得的结果(存储在内存中的数据)

_b_x1   _b_x2   _b_cons
.4470155    1.50748     1.043514
.4235979    1.60144     1.048863
.5006762    1.362679    .8828927
.5319981    1.494726    1.103693
.4926634    1.476443    .8611253
.5920001    1.557737    .8391003
.5893909    1.384571    1.312495
.4721891    1.37305     1.017576
.7109139    1.47294     1.055216
.4197589    1.442816    .9404677

上面Stata程序的R复制。

使用R，我设法得到以下信息(不是一个R专家)。然而，最让我担心的部分是循环结构，它遍历每一个重复的nreps的次数。

1. 定义reg_simulation函数。

#Defining a function 
reg_simulation<- function(obs = 1000){
    data <- data.frame(
    #Generate data
    x1 <-rnorm(obs, 0 , 1) ,
    x2 <-rnorm(obs, 0 , 1) ,
    y <- 1 + 0.5* x1 + 1.5 * x2 + rnorm(obs, 0 , 1) )
  #Estimate OLS
  ols <- lm(y ~ x1 + x2, data=data)  
  return(ols$coefficients)  
}

1. 使用for-循环结构调用reg_simulation 10次：

#Generate list to store results from simulation
results_list <- list()
# N repetitions
nreps <- 10
for (i in 1:nreps) {
  #Set seed internally (to get different values in each run)
  set.seed(i)
  #Save results into list
  results_list[i]  <- list(reg_simulation(obs=1000))  
}
#unlist results
df_results<- data.frame(t(sapply(results_list, 
                       function(x) x[1:max(lengths(results_list))])))

1. 结果：df_results。

#final results
df_results
#   X.Intercept.  x1        x2
# 1     1.0162384 0.5490488 1.522017
# 2     1.0663263 0.4989537 1.496758
# 3     0.9862365 0.5144083 1.462388
# 4     1.0137042 0.4767466 1.551139
# 5     0.9996164 0.5020535 1.489724
# 6     1.0351182 0.4372447 1.444495
# 7     0.9975050 0.4809259 1.525741
# 8     1.0286192 0.5253288 1.491966
# 9     1.0107962 0.4659812 1.505793
# 10    0.9765663 0.5317318 1.501162

Answer 1

你走在正确的轨道上。几个提示/更正：

1. 不要在<- data.frame()中使用

在R中，我们使用=构造内部列分配的数据帧，即data.frame(x = 1:10, y = 11:20)而不是data.frame(x <- 1:10, y <- 11:20)。

(有关于更多的话要说 <- vs =的报道，但我不想分散你的注意力。)

在您的示例中，您甚至不需要创建数据框架，因为x1、x2和y都将被识别为函数范围内的“全局”变量。我将在回答的末尾发布一些代码来演示这一点。

1. 通过R中的 for 循环生成列表时，始终尝试预先分配列表第一。

如果要增长(长) for循环，请始终尝试预先分配列表长度和类型。原因:这样，R知道有多少内存可以有效地分配给你的对象。在您只执行10次代表的情况下，这将意味着从以下内容开始：

results_list <- vector("list", 10)

lapply 3.考虑使用而不是 for

for循环在R中有一些不好的代表(这有点不公平，但这是另一天的故事)。许多R用户会考虑的另一种方法是lapply提供的函数式编程方法。我将暂停向您展示代码，但它看起来非常类似于for循环。请快速注意，从第2点开始，一个直接的好处是您不需要预先使用lapply.分配列表。

4.在并行中运行大循环

蒙特卡罗模拟是并行运行一切的理想选择，因为每次迭代都应该独立于其他迭代。在R中实现并行的一个简单方法是通过future.apply包。

把所有的东西放在一起，我可能会这样做你的模拟。请注意，这可能比你可能需要的更“先进”，但既然我在这里.

library(data.table)   ## optional, but what I'll use to coerce the list into a DT
library(future.apply) ## for parallel stuff
plan(multisession)    ## use all available cores

obs <- 1e3

# Defining a function 
reg_simulation <- function(...){
    x1 <- rnorm(obs, 0 , 1)
    x2 <- rnorm(obs, 0 , 1)
    y <- 1 + 0.5* x1 + 1.5 * x2 + rnorm(obs, 0 , 1)
    #Estimate OLS
    ols <- lm(y ~ x1 + x2)  
    
    # return(ols$coefficients)
    return(as.data.frame(t(ols$coefficients)))
}

# N repetitions
nreps <- 10

## Serial version
# results  <- lapply(1:nreps, reg_simulation)

## Parallel version
results  <- future_lapply(1:nreps, reg_simulation, future.seed = 1234L)

## Unlist / convert into a data.table
results <- rbindlist(results)

Answer 2

因此，在注释之后，您希望更改自变量(x)和误差项，并模拟系数，但如果出现错误，您也希望捕获错误。以下是其中的诀窍：

set.seed(42)
#Defining a function 
reg_simulation<- function(obs = 1000){

    data <- data.frame(
    #Generate data
    x1 <-rnorm(obs, 0 , 1) ,
    x2 <-rnorm(obs, 0 , 1) ,
    y <- 1 + 0.5* x1 + 1.5 * x2 + rnorm(obs, 0 , 1) )
  #Estimate OLS

    tryCatch(
      {
        ols <- lm(y ~ x1 + x2, data=data)  
        return(ols$coefficients)      
      }, 
      error = function(e){
              return(c('(Intercept)'=NA, 'x1'=NA, 'x2'=NA))
      }
    )
    
}
output <- t(data.frame(replicate(10, reg_simulation())))
output

    (Intercept)        x1       x2
X1    0.9961328 0.4782010 1.481712
X2    1.0234698 0.4801982 1.556393
X3    1.0336289 0.5239380 1.435468
X4    0.9796523 0.5095907 1.493548
...

在这里，tryCatch (也请参阅failwith)捕获错误并返回NA作为默认值。

请注意，您只需要设置一次种子，因为每次调用随机数生成器时，种子都会自动更改。