0

我正在为正态分布做一个循环，偏斜正态，而且偏斜很重。但是，我只能设法使用相同的 set.seed 获得正常和偏斜正常的 I 型错误。当我对 skewed heavy 运行编码时，我得到 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 值。

4

1 回答 1

0

我想你可能想在代码审查中发帖，但我最大的建议是更多地依赖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))
  }
)
于 2020-11-27T12:30:55.600 回答