測試哪個分佈有“長尾”
我測量了兩個非負隨機變量 A 和 B。它們的真實潛在概率是未知的，但是，可以假設概率在零時最大，並且對於較大的值單調遞減。最肯定的是，這些值也有一個上限，我對此有一個猜測，但這個猜測不是很好。
我想測試一下“尾巴”是否 $ P(A) $ 比尾巴“走得更遠” $ P(B) $ . 看起來確實如此，但也許這是偶然的？我可以考慮哪些指標？我試圖檢查平均值，但這兩個變量似乎具有可比性。
這個問題的基本特徵是：
- 它沒有做出強有力的分佈假設，因此具有非參數的味道。
- 它只涉及尾部行為，而不涉及整個分佈。
帶著一些信心——因為我沒有從理論上研究我的建議以完全理解它的性能——我將概述一種可能可行的方法。它藉鑑了 Kolmogorov-Smirnov 檢驗背後的概念、熟悉的基於秩的非參數檢驗和探索性數據分析方法。
讓我們從可視化問題開始。 我們可以在公共軸上繪製數據集的經驗分佈函數來比較它們：
黑色曲線顯示數據集 $ A $ （這裡有 $ m=50 $ 值），紅色曲線顯示數據集 $ B $ （這裡有 $ n=100 $ 值）。曲線在某個值處的高度 $ x $ 顯示值小於或等於的數據集的比例 $ x. $
這是一種情況，其中上半部分的數據 $ A $ 一直超過上半部分的數據 $ B. $ 我們可以看到，因為從左到右（從低值到高值）掃描，曲線最後交叉的高度為 $ 0.5 $ 之後，曲線為 $ A $ （黑色）保持在曲線的右側——即，高於——曲線的 $ B $ （紅色的）。這是數據分佈中右尾較重的證據 $ A $ 繪製。
**我們需要一個檢驗統計量。**它必須是一種以某種方式量化是否以及多少的方法 $ A $ 有一個“更重的右尾巴”比 $ B. $ 我的建議是這樣的：
- 將兩個數據集組合成一個數據集 $ n+m $ 價值觀。
- 對它們進行排名：這會分配值 $ n+m $ 到最高， $ n+m-1 $ 到下一個最高值，依此類推到該值 $ 1 $ 為最低。
- 加權排名如下：
- 劃分等級為 $ A $ 經過 $ m $ 和排名 $ B $ 經過 $ n. $
- 否定結果 $ B. $
- 累積這些值（以累積總和形式），從最大排名開始向下移動。
- 或者，通過將其所有值乘以某個常數來標準化累積和。
使用等級（而不是常數值 $ 1, $ 這是另一種選擇）對我們想要關注的最高值進行加權。該算法創建一個運行總和，當一個值來自 $ A $ 出現並且（由於否定）當值來自 $ B $ 出現。如果他們的尾巴沒有真正的區別，這個隨機遊走應該在零附近上下反彈。（這是加權的結果 $ 1/m $ 和 $ 1/n. $ ) 如果其中一條尾巴較重，則隨機遊走最初應向上趨向於較重 $ A $ 尾巴，否則頭向下以獲得較重 $ B $ 尾巴。
這提供了一個很好的診斷圖。 在圖中，我通過將所有值乘以來歸一化累積和 $ 1/\sqrt{n+m+1} $ 並按數字索引它們 $ q = 0/(m+n), 1/(m+n), \ldots, (m+n-1)/(m+n). $ 我稱之為“cranksum”（累積秩和）。這裡是前半部分，對應所有數據的上半部分：
有明顯的上升趨勢，與我們在上圖中看到的一致。 但這很重要嗎？
在零假設（同樣重尾）下的曲柄和模擬將解決這個問題。這樣的模擬創建了許多與原始數據相同大小的數據集 $ A $ 和 $ B $ （或者，幾乎等效地，創建組合數據集的許多任意排列）根據相同的分佈（哪個分佈無關緊要，只要它是連續的）；計算他們的曲柄和；並繪製它們。這是我為大小數據集製作的 40,000 個中的前 1000 個 $ 50 $ 和 $ 100: $
中間微弱的灰色鋸齒曲線形成一千個曲柄圖的組合。黃色區域，以粗體曲線（“信封”）為界，勾勒出鞋幫 $ 99.25 $ 和更低 $ 0.75 $ 所有 40,000 個值的百分位數。為什麼是這些百分位數？因為對這些模擬數據的一些分析表明*，在某些時候，*只有 5% 的模擬曲線超過了這些邊界。因此，因為實際數據的曲柄和圖確實超過了一些初始（低）值的上限 $ q, $ 它構成了重要的證據 $ \alpha=0.05 $ 水平（1）尾巴不同和（2）尾巴 $ A $ 比尾巴還重 $ B. $
**當然，您可以在圖中看到更多：**對於所有值，我們數據的曲柄和都非常高 $ q $ 之間 $ 0 $ 和 $ 0.23, $ 大約，然後才開始下降，最終達到 $ 0 $ 大約 $ q=0.5. $ 由此可見，至少上 $ 23% $ 數據集的底層分佈 $ A $ 持續超過上限 $ 23% $ 數據集的基礎分佈 $ B $ 並且可能是上層 $ 50% $ 的 … $ A $ 超過上限 $ 50% $ 的 … $ B. $
（因為這些是合成數據，我知道它們的基本分佈，所以我可以計算出對於這個例子，CDF 交叉在 $ x=1.2149 $ 在高度 $ 0.6515, $ 暗示上 $ 34.85% $ 的分佈為 $ A $ 超過 $ B, $ 非常符合曲柄和分析根據樣本告訴我們的內容。）
顯然，計算曲柄和並運行模擬需要一些工作，但它可以有效地完成：例如，這個模擬需要兩秒鐘。為了讓你開始，我已經附加了
R用於製作數字的代碼。# # Testing whether one tail is longer than another. # The return value is the cranksum, a vector of length m+n. # cranksum <- function(x, y) { m <- length(x) n <- length(y) i <- order(c(x,y)) scores <- c(rep(1/m, m), rep(-1/n, n)) * rank(c(x,y)) cumsum(scores[rev(i)]) / sqrt(n + m + 1) } # # Create two datasets from two different distributions with the same means. # mu <- 0 # Logmean of `x` sigma <- 1/2 # Log sd of `x` k <- 20 # Gamma parameter of `y` set.seed(17) y <- rgamma(100, k, k/exp(mu + sigma^2/2)) # Gamma data x <- exp(rnorm(50, mu, sigma)) # Lognormal data. # # Plot their ECDFs. # plot(ecdf(c(x,y)), cex=0, col="00000000", main="Empirical CDFs") e.x <- ecdf(x) curve(e.x(x), add=TRUE, lwd=2, n=1001) e.y <- ecdf(y) curve(e.y(x), add=TRUE, col="Red", lwd=2, n=1001) # # Simulate the null distribution (assuming no ties). # Each simulated cranksum is in a column. # system.time(sim <- replicate(4e4, cranksum(runif(length(x)), runif(length(y))))) # # This alpha was found by trial and error, but that needs to be done only # once for any given pair of dataset sizes. # alpha <- 0.0075 tl <- apply(sim, 1, quantile, probs=c(alpha/2, 1-alpha/2)) # Cranksum envelope # # Compute the chances of exceeding the upper envelope or falling beneath the lower. # p.upper <- mean(apply(sim > tl[2,], 2, max)) p.lower <- mean(apply(sim < tl[1,], 2, max)) # # Include the data with the simulation for the purpose of plotting everything together. # sim <- cbind(cranksum(x, y), sim) # # Plot. # q <- seq(0, 1, length.out=dim(sim)[1]) # The plot region: plot(0:1/2, range(sim), type="n", xlab = "q", ylab = "Value", main="Cranksum Plot") # The region between the envelopes: polygon(c(q, rev(q)), c(tl[1,], rev(tl[2,])), border="Black", lwd=2, col="#f8f8e8") # The cranksum curves themselves: invisible(apply(sim[, seq.int(min(dim(sim)[2], 1e3))], 2, function(y) lines(q, y, col="#00000004"))) # The cranksum for the data: lines(q, sim[,1], col="#e01010", lwd=2) # A reference axis at y=0: abline(h=0, col="White")
