生成不同变异性样本时检验统计量D1值相同的问题排查

问题描述 #

• 目标：生成两组不同变异性的样本，计算检验统计量D1，预期得到不同结果
• 问题：将new2_samples的生成项从sqrt(200)*new1_samples替换为sqrt(800)*new1_samples以提升变异性后，D1值完全相同
• 复现代码：

library(MASS)
set.seed(123)
# generate first sample
sample1 = replicate(2, rnorm(10,0,1), simplify=FALSE)
set.seed(456)
sample2 = replicate(2, rnorm(10,0,1), simplify=FALSE)

set.seed(789)
sample3 = replicate(2, rnorm(10,0,1), simplify=FALSE)
set.seed(101)
sample4 = replicate(2, rnorm(10,0,1), simplify=FALSE)
# express the data in matrix format
for (i in 1:2) {  
  sample1R = matrix(sample1[[i]], nrow=10, ncol=1)
  sample2R = matrix(sample2[[i]], nrow=10, ncol=1)
  
  # get the first set of data  
  sum_matrix = sample1R + sample2R
  new1_samples <- sqrt(0.5)* sum_matrix
  
  # set the mean to 100 and add it to the data
  mean_R=100
  # use constant c = sqrt(200 ) to reflect variation  of the data
  new2_samples<- sqrt(200)*new1_samples + mean_R
  
  #generate  second data  set
  sample1T = matrix(sample3[[i]], nrow=10, ncol=1)
  sample2T = matrix(sample4[[i]], nrow=10, ncol=1)
  
  Tsum_matrix = sample1T + sample2T
  Tnew1_samples <- sqrt(0.5)* Tsum_matrix
  
  # I will not multiply this data set with any constant to make sure it does not have same variabilities as the first data set  
  mean_T=100
   Tnew2_samples<- Tnew1_samples + mean_T
  
  # Calculate the standard deviation for each of the two data set 
  sd_T <- sd(Tnew2_samples)
  sd_R <- sd(new2_samples)
  
  # up_new2_row1 <- new2_samples[1] 
  up_new2_row1 <- new2_samples[1] + 5*sd_R
  
  up_Tnew2_row1 <- Tnew2_samples[1]
  #combine the two data sets  
  up_firstrow= cbind(up_new2_row1[],up_Tnew2_row1[]  )
  #combine  both full  columns 
  up_matrix= cbind(new2_samples,Tnew2_samples )
  
  # Update the up_samples matrix
  up_samples <- rbind(up_firstrow[],up_matrix[-1, ]  )
  
  # calculating my test stat
  H= up_samples
  mean_data1= mean(H[,1] )
  mean_data2= mean(H[,2] )
  
  mean_matrix=matrix( c(mean_data1, mean_data2  ), nrow=1, ncol=2, byrow=TRUE)
  
  residu1= H[1,]- mean_matrix
  cov_matrix= cov(H)
  
  n=10  
  s= (n-1)*cov(H)
  
  s_inverse<- solve(s)
  
  residu1_prime= t(residu1)
  # the test stat
  D1= residu1 %*% s_inverse %*% residu1_prime
  
  T_1= ((n-2)*D1)/((( n-1)/n)-D1 )
  
  cat("Simulation", i, ":
")
  
  print("up_samples")
  print(up_samples)
  
  print(" D1")
  print( D1)
  print(" T_1")
  print( T_1)
}

原因分析 #

D1是标准化后的马氏距离类统计量，你对new2_samples的缩放操作，被后续对第一个观测值的调整和协方差矩阵的缩放完全抵消了：

1. 当你将new2_samples乘以k（比如从sqrt(200)到sqrt(800)，k变为原来的2倍），sd_R也会同步变为原来的k倍
2. 第一个观测值up_new2_row1 = new2_samples[1] +5*sd_R，相当于被缩放了k倍（k*original_val +5*k*original_sd = k*(original_val+5*original_sd)）
3. 协方差矩阵cov(H)中，第一列的方差变为k²倍，与第二列的协方差变为k倍
4. 残差residu1 = H[1,]-mean_matrix的第一元素也被缩放了k倍
5. 协方差矩阵的逆矩阵s_inverse中，对应第一行第一列的元素变为1/k²倍，相关项变为1/k倍
6. 最终计算D1 = residu1 %*% s_inverse %*% t(residu1)时，k、1/k²、k的缩放因子完全抵消，结果保持不变

解决方案 #

要让D1随样本变异性改变，需要打破这种缩放抵消的逻辑，以下两种可行方式：

方式1：固定第一个观测值的增量（不基于缩放后的标准差） #

将up_new2_row1 <- new2_samples[1] +5*sd_R改为固定增量，比如：

up_new2_row1 <- new2_samples[1] + 5

这样修改sqrt(200)到sqrt(800)时，残差的变化不会被标准差的缩放抵消，D1会产生差异。

方式2：改变两组数据的相对方差比例 #

同时调整两组数据的缩放倍数，让它们的方差比发生变化，比如：

# 第一组用sqrt(200)缩放
new2_samples<- sqrt(200)*new1_samples + mean_R
# 第二组用sqrt(100)缩放（之前是不缩放，现在主动缩放改变比例）
Tnew2_samples<- sqrt(100)*Tnew1_samples + mean_T

这种情况下，两组数据的相对变异性改变，D1会随之变化。

修改后代码示例（方式1） #

library(MASS)
set.seed(123)
sample1 = replicate(2, rnorm(10,0,1), simplify=FALSE)
set.seed(456)
sample2 = replicate(2, rnorm(10,0,1), simplify=FALSE)

set.seed(789)
sample3 = replicate(2, rnorm(10,0,1), simplify=FALSE)
set.seed(101)
sample4 = replicate(2, rnorm(10,0,1), simplify=FALSE)

for (i in 1:2) {  
  sample1R = matrix(sample1[[i]], nrow=10, ncol=1)
  sample2R = matrix(sample2[[i]], nrow=10, ncol=1)
  
  sum_matrix = sample1R + sample2R
  new1_samples <- sqrt(0.5)* sum_matrix
  
  mean_R=100
  # 修改为sqrt(800)提升变异性
  new2_samples<- sqrt(800)*new1_samples + mean_R
  
  sample1T = matrix(sample3[[i]], nrow=10, ncol=1)
  sample2T = matrix(sample4[[i]], nrow=10, ncol=1)
  
  Tsum_matrix = sample1T + sample2T
  Tnew1_samples <- sqrt(0.5)* Tsum_matrix
  
  mean_T=100
  Tnew2_samples<- Tnew1_samples + mean_T
  
  sd_T <- sd(Tnew2_samples)
  sd_R <- sd(new2_samples)
  
  # 固定增量，不再乘以sd_R
  up_new2_row1 <- new2_samples[1] + 5
  up_Tnew2_row1 <- Tnew2_samples[1]
  
  up_firstrow= cbind(up_new2_row1[],up_Tnew2_row1[]  )
  up_matrix= cbind(new2_samples,Tnew2_samples )
  
  up_samples <- rbind(up_firstrow[],up_matrix[-1, ]  )
  
  # 计算检验统计量
  H= up_samples
  mean_data1= mean(H[,1] )
  mean_data2= mean(H[,2] )
  
  mean_matrix=matrix( c(mean_data1, mean_data2  ), nrow=1, ncol=2, byrow=TRUE)
  
  residu1= H[1,]- mean_matrix
  cov_matrix= cov(H)
  
  n=10  
  s= (n-1)*cov(H)
  
  s_inverse<- solve(s)
  
  residu1_prime= t(residu1)
  D1= residu1 %*% s_inverse %*% residu1_prime
  
  T_1= ((n-2)*D1)/((( n-1)/n)-D1 )
  
  cat("Simulation", i, ":
")
  print("up_samples")
  print(up_samples)
  print(" D1")
  print( D1)
  print(" T_1")
  print( T_1)
}

内容的提问来源于stack exchange，提问作者Amal Ghania