探索 t Statistic的稳健性

探索t-statistic 的稳健性
假设有2组数据，x1,…,xm; y1,…,yn。 t-test的假设为两组数据的平均值无差异 即：

令X和Sx为为x数组的平均值和标准差，Y和Sy为y数组的平均值和标准差；则H0的判别式为

在H0假设下，参数T符合t-分布（m+n-2*度）的前提：
1）数组x和数组y都是独立随机样品，符合正态分布；
2）两个群体的标准差相同
若显著性水平设置为a，则当 |T|≥tn+m-2,a/2,时 H0假设不成立。
什么是显著性水平a？？
显著性水平是估计总体参数落在某一区间内，可能犯错误的概率，用α表示；假设检验是围绕对原假设内容的审定而展开的，显著性是对差异的程度而言的，程度不同说明引起变动的原因也有不同：一类是条件差异，一类是随机差异。它是在进行假设检验时事先确定一个可允许的作为判断界限的小概率标准。
tn+m-2，a/2 即tdf，a 代表什么？代表在一个*度为df的t分布中，（1-a）所在的分位数。

##构造一个函数来计算t 值
#模拟构造一个符合正态分布的x和y数组，均值为50，标准差为10，一共10个数值
x=rnorm(10,mean=50,sd=10)
y=rnorm(10,mean=50,sd=10)
#获得数组长度
m = length(x)
n = length(y)
#计算Sp和t.stat
sp=sqrt(((m-1)*sd(x)^2+(n-1)*sd(y)^2)/(m+n-2))
t.stat=(mean(x)-mean(y))/(sp*sqrt(1/m+1/n))

# 通过构建函数，可以将功能打包
tstatistic=function(x,y)
{
  m=length(x)
  n=length(y)
  sp=sqrt(((m-1)*sd(x)^2+(n-1)*sd(y)^2)/(m+n-2))
  t.stat=(mean(x)-mean(y))/(sp*sqrt(1/m+1/n))
  return(t.stat)
}
# 若将函数另存为*.R文件，则可以直接调用
source("tstatistic.R")

#实际操作如下
data.x = c(1,4,3,6,5)
data.y = c(4,4,7,6,10,1212)
T = tstatistic(data.x,data.y)
# 已求得T值，如何判断是否存在差异？查表？？
t = qt(0.95,9)#计算n=5，m=6的t分布中，显著性水平（a=0.1）所在分位数的t值，即查表所得，判断|T|是否大于t，若小于t,则H0成立为大概率时间，小概率（<0.05）的情况未发生/
abs(T)

## 若 两个群体不符合正态分布和变异度相当的假设；则显著性水平则依赖于
  # 事先提供的 a 值
  # 数据群体的形状，例如normal, skewed, heavy-tailed等
  # 用两个标准差测量的两个群体的分布
  # the sample sizes m and n，群体的大小

#在给定 α, shape, spreads, and sample sizes,若需要估算αT（the true significance level）
    # αT = P(|T| ≥ tn+m−2,α/2)

#其基本算法是：
  #1 从初始群体中模拟x1-xm；y1-ym。
  #2 计算T值
  #3 判断|T|是否超过临界值，H0是否成立
  # 重复1-3步骤N次， 
    #αT = (number of rejections of H0)/N
#R 代码如下，其中alpha是给定的显著性水平，m，n为样本大小，qt(p,df)为在df*度的t分布中p分位数所对应的t值。
#abs(t) > qt(1-alpha/2,m+n-2)
#实际演示
alpha=.1; m=10; n=10 # 设置 alpha, m, n
N=10000 # 设置模拟次数
n.reject=0 # 计数器，累计H0不成立的次数

for (i in 1:N)
{
  x=rnorm(m,mean=0,sd=1) # 模拟样本x，基于均值0，sd1群体
  y=rnorm(n,mean=0,sd=1) # 模拟样本y，基于均值0，sd1群体
  t.stat=tstatistic(x,y) # 计算t值
  if (abs(t.stat)>qt(1-alpha/2,n+m-2))
    n.reject=n.reject+1 # 若H0不成立，即|T|> pt,则累加
}
true.sig.level=n.reject/N # 计算实际aT