如何编写拒绝 R 中的 Z 分数的函数

Question

当 zscore>=qnorm（1-alpha/2） 进行 10 次模拟且样本量为 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)

谢谢！

Answer 1

在您的编辑之后，我相信您要做的是查看在sims试验中，从正态分布中获取的大小n样本计算的偏度有多少次会被alph水平的显著性检验拒绝，因为偏差太大。

您有几个编码问题

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

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

   0.003623371

0.000000000 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")
        }
    }
}

但是，这个 stil 存在一些问题：

1. 从编程的角度来看
   • 打印出"拒绝"会立即提供反馈，但可扩展性不强。如果您有sims=1000，最好返回拒绝次数，nr 。如果您仍然想打印出"拒绝"nr次，您可以这样做:)
   • 此外，代码可以更简单，并且以更多的R风格编写，矢量化而不是使用循环。这也将具有更快的优点。由于 R 是一种解释型语言，因此矢量化会产生巨大的差异，因为数字处理可以在后台进行，而无需 R 一遍又一遍地遍历for循环。
2. 也许更严重的是，存在一些统计问题：
   • 6/n 是偏度方差的估计值（维基百科），但您需要标准差，因此需要取 6/n 的平方根。
   如果 z 分数大于第 1-alph/2 个分位数，则代码拒绝，
   • 但如果 z 分数小于第 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)
}

您应该能够修改它以适应您的目的，但如有必要，我可以澄清。

Answer 2

不确定您要实现什么，但这里有一种方法可以做到这一点

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 

Answer 3

您错误地循环了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)