显示 1000 次模拟的 β̂ 1 OLS 和 β̂ 1 ROLS 的核密度估计

Posted 2023-02-16

【中文标题】显示 1000 次模拟的 β̂ 1 OLS 和 β̂ 1 ROLS 的核密度估计

【英文标题】：showing the kernel density estimates of β̂ 1 OLS and β̂ 1 ROLS for your 1000 simulations

【发布时间】：2022-01-21 02:20:37

【问题描述】：

我有以下问题：我做了 1000 次模拟以获得 ß^1_OLS 和 ß^1_ROLS 的向量，现在我必须通过将它们都绘制为密度函数来比较这两者（我必须使用 ggplot( ））。是否可以根据向量绘制密度？

按照我的代码，请原谅混乱：

>N=10000
X <- runif(N, min = 0, max = 100)
E <- runif(N, min = -5, max = 5)

U = E * sqrt(X)

# population regression
Y <- 3 + 2 * X + U
population <- data.frame(X, Y)


# set sample size an repetition 
n <- 1000
reps<- 1000



# initialize the matrix of outcomes
fit_coef <- matrix(ncol = 2, nrow = reps)
fit_coef_ROLS<-matrix(ncol = 2, nrow = reps)


# #######loop sampling and estimation of the coefficients for OLS
set.seed(1)
for (i in 1:reps)
  
  sample <- population[sample(1:N, n), ]
  fit_coef[i,] <- lm(Y ~ X, data = sample)$coefficients
  

fit_coef_ß1<- fit_coef[,-1] 
fit_coef_ß1
  
#######loop sampling and estimation of the robust coefficients ROLS
set.seed(1)
for (i in 1:reps)
  
  sample <- population[sample(1:N, n), ]
  fit_coef_ROLS[i,] <- rlm(Y ~ X, data = sample)$coefficients
  

fit_coef_ß1_ROLS<- fit_coef_ROLS[,-1] 
fit_coef_ß1_ROLS
 

## Plot
df_coef_OLS<-as.data.frame(fit_coef)

plot_coef_OLS<-ggplot(df_coef_OLS, aes(x = fit_coef_ß1)) +
  geom_density()
plot_coef_OLS

如果我使用这个公式，我会得到一个图，但我不知道如何在 1 个图中获得两个估计量的密度，而且密度不在 [0,1] 的范围内

感谢您的支持！ This is the task

【问题讨论】：

密度不必小于1。必须为正但可以任意大。

【参考方案1】：

当然，您只需要重塑您的数据以具有长格式，即每次迭代的相关系数在 data.frame 中获得一行。然后，您需要另一列来声明该行来自哪个模拟。

library(ggplot2)

## Your code goes here (omitted for brevity)

df <- data.frame(
  coef = c(fit_coef[, -1], fit_coef_ROLS[, -1]),
  type = rep(c("OLS", "ROLS"), c(nrow(fit_coef), nrow(fit_coef_ROLS)))
)

ggplot(df, aes(coef, colour = type)) +
  geom_density()

^{由reprex package (v2.0.1) 于 2021 年 12 月 18 日创建}

【讨论】：

非常感谢！密度不在[0,1]之间是正常的还是还是有错误？

这很常见。回想一下密度积分为 1，因此密度是多少取决于数据的范围。在这里，数据（大约）在 1.95 和 2.05 之间，因此它应该平均在 10 左右才能积分为 1。相反，如果范围是，例如，在 0 和 10 之间，那么密度应该平均为 0.1 左右。跨度>