0

当 zscore>=qnorm(1-alpha/2) 对于 alpha=0.05 和样本大小为 10 的 10 次模拟时，我如何编写一个返回“拒绝”的函数。我编写了以下代码，但我没有得到正确的输出。“zscore”是检验统计量，t 是正常的，平均值和标准差为 6/n。sims 对应于要执行的模拟次数。此函数应模拟蒙特卡洛评估。

testsk=function(n,alph,sims){
    t=numeric(sims)
for (i in 1:sims) {
    x=rnorm(n)
t[i]=skewness(x)
zscore=t/(6/n)
return(zscore)
}
if(zscore>=qnorm(1-alph/2)){
print("REJECT")
}
}


testsk(10,0.05,10)

谢谢！

4

3 回答 3

1

Not sure what are you trying to achieve but here is one way you can do this

testsk <- function(n, alph, sims){
  for (i in 1:sims){
  x <- rnorm(n)
  zscore <- skewness(x)/(6/n)
  cat(paste0("Simulation #", i,":"), ifelse(zscore >= qnorm(1 - alph/2), "REJECT", "Don't REJECT"), "\n")
  }
}

n <- 10
alph <- .05
sims <- 10
testsk(n, alph, sims)

#Simulation #1: Don't REJECT 
#Simulation #2: REJECT 
#Simulation #3: Don't REJECT 
#Simulation #4: Don't REJECT 
#Simulation #5: Don't REJECT 
#Simulation #6: Don't REJECT 
#Simulation #7: Don't REJECT 
#Simulation #8: Don't REJECT 
#Simulation #9: Don't REJECT 
#Simulation #10: Don't REJECT
于 2014-04-01T22:37:11.940 回答
1

在您的编辑之后，我相信您想要做的是查看从正态分布中提取sims的大小样本计算出的偏度经过多少次试验，因为该水平的显着性检验过于偏斜而被拒绝。nalph

你有几个编码问题

  1. 您想为每次试验进行 z 分数测试，但测试在循环之外。
  2. z 分数是使用向量计算的t，但您想使用标量计算它t[i]。

  3. 循环内有一条return语句将导致函数在循环的第一次迭代中终止，并返回 z 分数。由于第 2 个原因，z 分数是一个向量，但它的倒数第二个值全为零，因为您只运行了一次迭代，因此该函数的典型输出如下

    0.003623371 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000

修复这些直接问题会产生以下代码

library(e1071)
testsk=function(n,alph,sims) {
    t=numeric(sims)
    for (i in 1:sims) {
        x=rnorm(n)
        t[i]=skewness(x)
        zscore=t[i]/(6/n)
        if(zscore>=qnorm(1-alph/2)){
            print("REJECT")
        }
    }
}

然而，这仍然存在一些问题：

  1. 从编程的角度来看
    • 打印出“REJECT”会立即提供反馈，但可扩展性不强。如果你有sims=1000你最好返回拒绝的数量，nr. 如果您仍然想打印“拒绝”nr时间，您可以这样做:)
    • 此外，代码可以更简单，并以更多的 R 风格编写，矢量化而不是使用循环。这也将具有更快的优势。因为 R 是一种解释性语言，所以向量化会产生巨大的差异，因为数字运算可以在引擎盖下发生，而 R 不必for一遍又一遍地遍历你的循环。
  2. 也许更严重的是，存在一些统计问题：
    • 6/n 是对偏度方差的估计（维基百科），但您需要标准差，因此您需要取 6/n 的平方根。
    • 如果 z 分数大于第 th 分位数，代码将拒绝，但如果 z 分数小于第1-alph/2th 分位数，它也应该拒绝alph/2。就目前而言，您的拒绝区域alph/2不是alph。
    • 可能还有其他问题，但我认为这些是主要问题。（我假设您知道 6/n 只是对大样本方差的一个很好的估计。）

正确的程序如下

library(e1071)
testsk=function(n,alph,sims) {
    # Generate random numbers in a matrix, each trial is a row
    X=matrix(rnorm(sims*n), ncol=n)
    # Get skewnesses, 1 means apply to rows
    skews=apply(X,1,skewness)
    # Calculate z score vector and rejection vector
    zscore=skews/sqrt(6/n)
    reject=!(qnorm(alph/2) < zscore & zscore < qnorm(1-alph/2))
    # Return the number of rejections
    sum(reject)
}

您应该能够对其进行修改以适合您的目的，但如有必要，我可以澄清一下。

于 2014-04-01T23:42:28.930 回答
0

您错误地循环了sims. 你能解释一下你想要做什么吗？

testsk <- function(n,alph,sims) {
  t <- numeric(sims)
  for (i in seq_along(sims)) {
    x <- rnorm(n)
    t[i] <- skewness(x)
  }
  zscore <- t/(6/n)
  if (any(zscore>=qnorm(1-alph/2))) {
    print("REJECT")
  }
  return(zscore)
}

testsk(10,0.05,10)
于 2014-04-01T22:03:03.873 回答