如何編寫一個拒絕R中Z分數的函數？

Question

如何寫一個函數返回「拒絕」當zscore> = qnorm（1-alpha/2）10個模擬alpha = 0.05和樣本大小爲10.我寫了下面的代碼，但我沒有獲取相關輸出。 「zscore」是檢驗統計量，t是平均值和標準偏差6/n的正態分佈。 sims對應於要執行的模擬的數量。該函數應該模仿Monte Carlo評估。

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試驗根據從正態分佈獲得的樣本計算的偏度將被拒絕，因爲在alph水平上的顯着性檢驗也偏斜。

你有幾個編碼問題

1. 您想爲每一個審判做Z值測試，但測試是外循環。
2. z分數是使用向量t計算得出的，但是您想使用標量t[i]來計算z分數。

3. 循環內部有一個return語句，它將導致函數在循環的第一次迭代中終止，返回z分數。對於沒有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是偏度（wikipedia）的方差的估計，但你想要的標準偏差，所以你需要採取的平方根6/n。
   • 如果z分數大於1-alph/2 th分位數，則代碼拒絕，但如果z分數小於alph/2 th分位數，它也應拒絕。因爲它代表你的拒絕區域是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

您錯誤地循環了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) 

Answer 3

不知道什麼是你想實現，但在這裏是一種方法可以做到這一點

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