如何解决使用 R 进行蒙特卡洛模拟如何为 R“翻译”Stata 代码
我需要从 Stata 到 R 运行相同的练习(蒙特卡罗模拟)。
我在 Stata 中使用的代码是下面的代码。我如何使用 R 做到这一点? (我已经搜索了很多教程,但我仍然没有设法在 R 中做到这一点)。
* Simulations (10,100 and 1000 sample replications/iterations)
clear
drop _all
set obs 100
set seed 54231
gen x = ((rnormal()))*10 + 40
* Generating true_y,considering Beta = 0,035
gen true_y = 5+0.03500*x
save truth,replace
twoway scatter true_y x
program hetero1
version 13
args c
use truth,clear
gen y = true_y + rnormal()
regress y x
end
foreach i in 10 100 1000 {
simulate _b,reps (`i'): hetero1
sum _b_x
twoway histogram _b_x,fraction xline(+0.03500,lcolor(red)) xline(`r(mean)',lcolor(green)) fcolor(none) lcolor(gs0) legend(off) title(`i' Repetições)
graph save graf`i'.gph,replace
}
gr combine graf10.gph graf100.gph graf1000.gph
graph export "graf.png",as(png) replace
最后,考虑到 10、100 和 1000 个样本重复,我需要获得这 3 个直方图(估计的 beta/系数)。红线指的是“真实”系数,绿线是估计系数的平均值 - [见链接中的图片]
解决方法
应该这样做:
# Simulations (10,100 and 1000 sample replications/iterations)
library(ggplot2)
library(dplyr)
library(gridExtra)
n <- 100
set.seed(54231)
x <- rnorm(n)*10 + 40
# Generating true_y,considering Beta = 0,035
true_y <- 5+0.03500*x
plot(x,true_y)
b <- t(replicate(1110,coef(lm(true_y + rnorm(n) ~ x))))
b <- as.data.frame(b) %>%
rename("a" = "(Intercept)","b" = "x") %>%
mutate(
obs = 1:n(),n = case_when(
obs %in% 1:10 ~ "N = 10",obs %in% 11:110 ~ "N = 100",TRUE ~ "N = 1000"),n = factor(n,levels=c("N = 10","N = 100","N = 1000")))
b10 <- b %>% filter(n == "N = 10")
g1 <- ggplot() +
geom_histogram(data = b10,aes(x=b),bins=3,col="white") +
geom_vline(xintercept = 0.03500,col="red") +
geom_vline(data = b10 %>% summarise(b=mean(b)),aes(xintercept = b),col="green") +
facet_wrap(~n) +
theme_bw()
b100 <- b %>% filter(n == "N = 100")
g2 <- ggplot() +
geom_histogram(data = b100,bins=10,col="red") +
geom_vline(data = b100 %>% summarise(b=mean(b)),col="green") +
facet_wrap(~n) +
theme_bw()
b1000 <- b %>% filter(n == "N = 1000")
g3 <- ggplot() +
geom_histogram(data = b1000,bins=25,col="red") +
geom_vline(data = b1000 %>% summarise(b=mean(b)),col="green") +
facet_wrap(~n) +
theme_bw()
library(gridExtra)
grid.arrange(g1,g2,g3,nrow=2)
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。