R上的随机样本和分布

bzzcjhmw 于 2022-12-06 发布在其他

关注(0)|答案(2)|浏览(148)

我试图建立一个估计量来比较GLS估计量和OLS估计量的渐近性，我的想法是尝试看看在大样本情况下会发生什么。
理想情况下，我希望创建一个循环，该循环将为不同的参数值生成6000个不同的随机样本，每个样本的大小为50和100。

N=1000
n=c(50, 100)

#parameters
alpha0=1
beta0=1
gamma0=c(0, 0.1, 0.5)

alpha1=matrix(NA,N,6)
beta1=matrix(NA,N,6)
alpha2=matrix(NA,N,6)
beta2=matrix(NA,N,6)
alphaOLS=matrix(NA,N,6)
betaOLS=matrix(NA,N,6)

不同的样本来自gamma 0和n的组合，这将等于6（乘以N），得到6000。我的第一个想法是构建一个用于生成随机样本的循环，我尝试使用的模型如下：y_i=alpha+betax_i+u_i
u_i=e_ih（x_i）^（1/2）且h（x）=exp（γ 0）

u <- list()

for (i in n) {
  for (k in gamma0) {
    x=rnorm(i,0,1)
    h=exp(gamma0[k]*x)
    e=rnorm(i,0,1)
    u[[i]] <- e*h^(1/2)
    
  }
  
}

这个循环的问题是，我只得到x和e中的一个随机样本，h作为空矩阵出来，因此，u也是空的。h这里应该是一个矩阵，其中列对应于x* γ 0的不同值。e应该是N（0，1），u是模型的残差
我的理想输出应该是让这个循环工作，因为从那里开始，我可以围绕手动构建OLS和GLS估计器进行排序。
多谢了！

来源：https://stackoverflow.com/questions/74574895/random-samples-and-distributions-in-r

2条答案

按热度按时间

tmb3ates1#

这应该可以工作，这里我们直接使用i和k，而不是循环中的n和gamma0。

### parameters
N <- 1000
n <- c(50, 100)

alpha0 <- 1
beta0 <- 1
gamma0 <- c(0, 0.1, 0.5)

### Initiating objects for the loop
u <- list()
num_iter <- 0

### Looping
for (i in n) {
  for (k in gamma0) {
    
    num_iter <- num_iter + 1
    
    x <- rnorm(i, 0, 1)
    h <- exp(k * x)
    
    e <- rnorm(i, 0, 1)
    u[[num_iter]] <- e * h^(1/2)
    names(u)[num_iter] <- paste("n:", i, ", gamma:", k, sep="", collapse=" ")
    
  }
}

### Display results
u

编辑（根据您最近的请求）

### Parameters for each scenario
N <- 1000
n <- c(50, 100)
gamma0 <- c(0, 0.1, 0.5)

### Initiating dataframes for each scenario
for (i in n) {
  for (k in gamma0) {
assign(paste("df_n", i, "gamma", k, sep=""), list())
  }
}

### Defining a dataframe with n rows (sample size) and N columns(simulation) for each of the 6 scenarios
df_n100gamma0 <- sapply(1:1000, function(x) rnorm(100, 0, 1) * exp(0 * rnorm(100, 0, 1))^(1/2))
df_n100gamma0.1 <- sapply(1:1000, function(x) rnorm(100, 0, 1) * exp(0.1 * rnorm(100, 0, 1))^(1/2))
df_n100gamma0.5 <- sapply(1:1000, function(x) rnorm(100, 0, 1) * exp(0.5 * rnorm(100, 0, 1))^(1/2))
df_n50gamma0 <- sapply(1:1000, function(x) rnorm(50, 0, 1) * exp(0 * rnorm(50, 0, 1))^(1/2))
df_n50gamma0.1 <- sapply(1:1000, function(x) rnorm(50, 0, 1) * exp(0.1 * rnorm(50, 0, 1))^(1/2))
df_n50gamma0.5 <- sapply(1:1000, function(x) rnorm(50, 0, 1) * exp(0.5 * rnorm(50, 0, 1))^(1/2))

### Dimension of 1 dataframe
dim(df_n100gamma0)

### Results
# First simulation
df_n100gamma0[, 1]

# 57th simulation
df_n100gamma0[, 57]

赞(0）回复(0）举报 2022-12-06

1szpjjfi2#

或者，通过使用R固有的矢量化和tidyverse处理组的能力来完全避免循环。使用R时的一个很好的经验法则是“如果我想使用for循环，可能有更好的方法来做到这一点”。
这里，我只是为了简单起见而使用OLS，但是使用这两种方法的扩展应该是直接的。您可以将代码 Package 在函数中，并为N、n和gamma0提供参数，以提供更大的通用性。

library(tidyverse)
library(broom)

# Create a tibble containing all combinations of the various design parameters
tibble() %>% 
       expand(
         # Sample=1:1000,
         # For speed
         Sample=1:10,
         Gamma0=c(0, 0.1, 0.5),
         n=c(50, 100),
         alpha0=1,
         beta0=1
       ) %>% 
       #  For each combination of design parameters...
       group_by(Gamma0, n, Sample, alpha0, beta0) %>% 
       group_map(
         function(.x, .y) {
           #  ... create some data.
           df <- .x %>% 
                    expand(
                      nesting(Gamma0, n, Sample, alpha0, beta0), 
                      ID=1:n
                    ) %>% 
                    mutate(
                      X=rnorm(n),
                      H=Gamma0 * X,
                      E=rnorm(n),
                      U=E*H^(1/2),
                      Y=alpha0 + beta0 * X + U
                    )
           # And analyse it. For example
           glm(Y ~ X, data=df) %>% 
             tidy() %>% 
             bind_cols(.y) %>% 
             add_column(model="OLS")
         },
         .keep=TRUE
       ) %>% 
       # Combine all analyses of each individual sample into a single tibble
       bind_rows()
# A tibble: 120 × 11
   term        estimate std.error statistic p.value Gamma0     n Sample alpha0 beta0 model
   <chr>          <dbl>     <dbl>     <dbl>   <dbl>  <dbl> <dbl>  <int>  <dbl> <dbl> <chr>
 1 (Intercept)     1     7.39e-17   1.35e16       0      0    50      1      1     1 OLS  
 2 X               1     6.79e-17   1.47e16       0      0    50      1      1     1 OLS  
 3 (Intercept)     1     0        Inf             0      0    50      2      1     1 OLS  
 4 X               1     0        Inf             0      0    50      2      1     1 OLS  
 5 (Intercept)     1.00  8.55e-17   1.17e16       0      0    50      3      1     1 OLS  
 6 X               1     8.37e-17   1.19e16       0      0    50      3      1     1 OLS  
 7 (Intercept)     1     4.00e-17   2.50e16       0      0    50      4      1     1 OLS  
 8 X               1     3.56e-17   2.81e16       0      0    50      4      1     1 OLS  
 9 (Intercept)     1     4.22e-17   2.37e16       0      0    50      5      1     1 OLS  
10 X               1     4.40e-17   2.27e16       0      0    50      5      1     1 OLS  
# … with 110 more rows

赞(0）回复(0）举报 2022-12-06

我来回答

R上的随机样本和分布

2条答案

相关问题

热门标签

最新问答