时间戳

首页 > 解决方案 > 等价于 R 中用于蒙特卡罗模拟的 Stata 命令“模拟”

r - 等价于 R 中用于蒙特卡罗模拟的 Stata 命令“模拟”

问题描述

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

下面是该simulate程序用法的 Stata 说明，以及我尝试使用R.

最后，我的主要问题是：R 用户将如何运行蒙特卡罗模拟？还是我在结构或速度瓶颈方面遗漏了什么？非常感谢您。

统计示例

  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复制。

使用 RI 已成功获得以下结果（不是 R 专家）。然而，最让我担心的部分是循环遍历每个重复次数的 for 循环结构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. reg_simulation使用 for 循环结构调用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

标签： rrandomsimulationstatamontecarlo

解决方案

因此，根据评论，您希望改变自变量 (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 作为默认值。

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

推荐阅读