在 R 中为正态分布设置种子

Question

我正在为正态分布、偏正态分布和偏重分布做循环。但是，我只能设法使用相同的 set.seed 获得正常和偏斜正常的 I 型错误。当我运行偏重的编码时，我得到NA，并且在编码中计算的数据中也有一些缺失值。我怎样才能解决这个问题？无论如何，我已经尝试了很多值而不是 val，但它并没有让我对正常和偏斜正常感到满意，因为我总是得到零。

asim<-10000
pv<- rep(NA,asim)

for (val in 1:asim){
  
  set.seed(val)
  
  t=3
  n1=10
  n2=10
  n3=10
  N=n1+n2+n3
  
  data1<-rnorm(n1,0,1)
  data2<-rnorm(n2,0,1)
  data3<-rnorm(n3,0,1)
  
  # SKEWED HEAVY
  
  G=0.5
  h=0.5 # h=0 if skewed normal
  
  y1=((exp(G*data1)-1)/G)*(exp((h*data1^2)/2))
  y2=((exp(G*data2)-1)/G)*(exp((h*data2^2)/2))
  y3=((exp(G*data3)-1)/G)*(exp((h*data3^2)/2))
  
       g1=sort(y1)
  g2=sort(y2)
  g3=sort(y3)
  
  
  ybar1<-mean(g1)
  ybar2<-mean(g2)
  ybar3<-mean(g3)
  
  
  #BIWEIGHT
  
  med1=median(g1)
  med2=median(g2)
  med3=median(g3)
  
  mad=mad(c(g1,g2,g3))
  
  u1=(g1-med1)/(9*mad)
  u2=(g2-med2)/(9*mad)
  u3=(g3-med3)/(9*mad)
  
  idx1<- which(abs(u1)<1)
  idx2<- which(abs(u2)<1)
  idx3<- which(abs(u3)<1)
  
  
  #create empty list
  
  num1=c()
  num2=c()
  num3=c()
  den1=c()
  den2=c()
  den3=c()
  nume1=c()
  nume2=nume3=deno1=deno2=deno3=c()
  
  
  for(z in idx1){
    num1[z] = ((g1[z]-med1)^2)*((1-(u1[z]^2))^4)
    den1[z] = ((1-(u1[z]^2))*(1-5*(u1[z]^2)))
  }
  
  for(j in idx2){
    num2[j] = ((g2[j]-med2)^2)*((1-(u2[j]^2))^4)
    den2[j] = ((1-(u2[j]^2))*(1-5*(u2[j]^2)))
  }
  
  for(k in idx3){
    num3[k] = ((g3[k]-med3)^2)*((1-(u3[k]^2))^4)
    den3[k] = ((1-(u3[k]^2))*(1-5*(u3[k]^2)))
  }
  
  print(num1)

 }

打印num1后，有NA值

Answer 1

我认为您可能希望在代码审查中 post，但我的最大建议是更多地依赖 lapply 而不是按顺序命名的对象（例如 num1, num2, num3.

我假设在此结束时，您需要模拟的结果。为此，我们将在每次模拟中使用 lapply 制作一个列表。我在循环之外添加了一些常量；我们需要在每次模拟期间复制 3 个数据集。我们广泛使用 lapply 来不断制作您需要的汇总统计数据。 Map 类似，只是它允许我们在需要时传递多个参数。

asim<-10
t=3
n = rep(10, t)
N = sum(n)
pv<- rep(NA,asim)

lapply(seq_len(asim), function (val) {
  set.seed(val)
  
  data = lapply(n, rnorm, 0, 1)
  
  # SKEWED HEAVY
  
  G=0.5
  h=0.5 # h=0 if skewed normal
  
  y = lapply(data, function(x) ((exp(G*x)-1)/G)*(exp((h*x^2)/2)))  
  
  g = lapply(y, sort)
  ybar = lapply(y, mean)
  
  #BIWEIGHT
  
  med = lapply(g, median)
  mad_=mad(unlist(g, use.names = FALSE))
  u = Map(function (g_, med_) (g_ - med_) / (9 * mad_), g, med)
  idx = lapply(u, function(x) which(abs(x) < 1))


  num = Map(function(g_, med_, u_, idx_) (((g_ - med_) ^ 2) * ((1 - (u_^2))^4))[idx_], g, med, u, idx)
  den = Map(function(u_, idx_) ((1 - (u_^2)) * (1 - 5 * (u_ ^ 2)))[idx_], u, idx)

  data.frame(sim = rep(seq_len(t), lengths(idx)), 
             num = unlist(num, use.names = FALSE), 
             den = unlist(den, use.names = FALSE))
  }
)

在 R 中为正态分布设置种子

Set seed in R for Normal distribution

r

random-seed