多元正态分布的均值向量的检验及R实现

$p$维正态总体$N_P(\mu,\sum)$的均值向量检验，$X_1,X_2,\cdots,X_n$是来自正态总体的样本：

1.$\sum$已知时单个总体均值向量的检验：

具体步骤：

• 作统计假设：$H_0:\mu=\mu_0,H_1:\mu\neq \mu_0$
• 计算样本的均值
• 计算统计量的具体值:
  $$\begin{cases} U=\frac{\bar{X}-\mu_0}{\sigma}\sqrt{n} & \text{ if } p=1 \\ {T_0}^2= n(\bar{X}-\mu_0)'\sum^{-1}(\bar{X}-\mu_0)& \text{ if } p>1 \end{cases}$$
• 按规定的小概率标准$\alpha$，查$\chi^2$分布表，得临界值$\chi_\alpha^2(p)$，并做出判断：
  -当${T_0}^2\leqslant \chi_\alpha^2(p)$，接受原假设，即认为没有显著差异
  -当${T_0}^2> \chi_\alpha^2(p)$，拒绝原假设，即认为有显著差异

*R实现：

mu.test.known=function(data, mu0, Sigma0, alpha=0.05)   
###############################################################
## H0: mu=mu0 when Sigma0 is known
## this is a Chisq testing
##############  Input  ########################################
## data  = design matrix with the ith sample in the ith line
## mu0   = mu0 for null hypothesis
## Sigma0= the known variance matrix
## alpha = the significant level, default value = 0.05
############## Output  ########################################
## Reject.area = reject region
## p.value     = p value
###############################################################
{
  data=as.matrix(data) #将数据框转化为矩阵#
  n=nrow(data) #n行#
  p=ncol(data) #p列#

  X.bar=apply(data, 2, mean) #按列求均值#
  T1=n*t(X.bar-mu0)%*%solve(Sigma0)%*%(X.bar-mu0)

  a2=qchisq(1-alpha, p) #求下侧分位点，上侧：lower.tail = FALSE#

  reject=matrix(c(T1, a2), nrow=1) #按行排#
  rownames(reject)=c("Reject") #行名#
  colnames(reject)=c("Obs", ">1-alpha") #列名#

  pv=1-pchisq(T1, p) #右半累积概率，T越大，P越小，越拒绝#
  return(list(Reject.area=reject, p.value=pv))
}

2.$\sum$未知时单个总体均值向量的检验：

具体步骤：

• 作统计假设：$H_0:\mu=\mu_0,H_1:\mu\neq \mu_0$
• 计算样本的均值$\bar{X}$和样本协方差$V=\frac{1}{n-1}\sum_{i=1}^{n}(X_i-\bar{X})(X_i-\bar{X})'$
• 计算统计量的具体值:
  $$T^2=n(\bar{X}-\mu_0)'V^{-1}(\bar{X}-\mu_0),\frac{n-p}{(n-1)p}T^2\sim F(p,n-p)$$
• 按规定的小概率标准$\alpha$，查$F$分布表，得临界值$F_\alpha(p,n-p)$，并做出判断：
  -当${F_0}\leqslant F_\alpha(p,n-p)$，接受原假设，即认为没有显著差异
  -当${F_0}> F_\alpha(p,n-p)$，拒绝原假设，即认为有显著差异

*R实现：

mu.test=function(data, mu0)   
###############################################################
## H0: mu=mu0 when Sigma is unknown
## this is an F testing
##############  Input  ########################################
## data  = design matrix with the ith sample in the ith line
## mu0   = mu0 for null hypothesis
############## Output  ########################################
## p.value     = p value
###############################################################
{
  data=as.matrix(data)
  n=nrow(data)
  p=ncol(data)

  X.bar=apply(data, 2, mean)
  A=(n-1)*var(data)

  T2=(n-1)*n*t(X.bar-mu0)%*%solve(A)%*%(X.bar-mu0)
  F=(n-p)/((n-1)*p)*T2

  p.two=1-pf(F, p, n-p)
  return(list(p.value=p.two))
}

补充：R生成多元正态分布

######################生成多元正态分布####################
library(MASS)
Sigma <- matrix(c(10,3,3,2),2,2)
Sigma
mvrnorm(n=1000, rep(0, 2), Sigma) 
##########################模拟##########################
library(MASS)
source("Norm Mean Test.r")
error=0
alpha=0.05
Sigma0 = matrix(c(1,0.5,0.5,1),2,2)
for(i in 1:1000){
data=mvrnorm(n=200, rep(0, 2), Sigma0) 
mu0=c(0,0)
result1=mu.test.known(data, mu0, Sigma0, alpha=0.05)
if (result1[2]<alpha) error=error+1
}
EERROR=error/1000

3.$\sum$未知但相等时两个总体均值向量的检验：

具体步骤：

• 作统计假设：$H_0:\mu_1=\mu_2,H_1:\mu_1\neq \mu_2$
• 计算样本的均值$\bar{X}$和$\bar{Y}$，样本离差阵$A_1=\sum_{i=1}^{n}(X_i-\bar{X})(X_i-\bar{X})'$和$A_2=\sum_{i=1}^{n}(Y_i-\bar{Y})(Y_i-\bar{Y})'$
• 计算统计量的具体值:
  $$T^2=\frac{nm}{n+m}(\bar{X}-\bar{Y})'V_e^{-1}(\bar{X}-\bar{Y})$$ 其中$V_e=\frac{1}{n+m-2}(A_1+A_2)$ $$\frac{n+m-p-1}{(n+m-2)p}T^2\sim F(p,n+m-p-1)$$
• 按规定的小概率标准$\alpha$，查$F$分布表，得临界值$F_\alpha(p,n+m-p-1)$，并做出判断：
  -当${F_0}\leqslant F_\alpha(p,n+m-p-1)$，接受原假设，即认为没有显著差异
  -当${F_0}> F_\alpha(p,n+m-p-1)$，拒绝原假设，即认为有显著差异

*R实现：

#### two independent normal distribution  ##################
two.mu.test=function(data1, data2)   
###################################################################
## H0: mu1=mu2 
## this is an F testing
##############  Input  ############################################
## data1  = design matrix for X with the ith sample in the ith line
## data2  = design matrix for X with the ith sample in the ith line
############## Output  ####################################### p.value     = p value
############################################################
{
  data1=as.matrix(data1)
  data2=as.matrix(data2)
  n1=nrow(data1)
  n2=nrow(data2)
  p=ncol(data1)

  X.bar=apply(data1, 2, mean) 
  A1=(n1-1)*var(data1)
  Y.bar=apply(data2, 2, mean)
  A2=(n2-1)*var(data2) 
  A=(A1+A2)/(n1+n2-2)

  T2=(n1*n2/(n1+n2))*t(X.bar-Y.bar)%*%solve(A)%*%(X.bar-Y.bar)
  F=(n1+n2-2-p+1)/((n1+n2-2)*p)*T2

  p.two=1-pf(F, p, (n1+n2-p-1))
  return(list(p.value=p.two))
}

4.多个正态总体均值向量的检验-多元方差分析：

（$k$个$N_p(\mu_t,\sum)$，其中$t=1,\cdots,k$，样本$X_{(i)}^{(t)}(t=1,\cdots,k)(i=1,\cdots,n_i)$)
实验数据：

水平数	重复数	平均
$X_1(p\times n_1)$	$X_{(1)}^{(1)},\cdots,X_{(n_1)}^{(1)}$	$\bar{X}^{(1)}$
$X_2(p\times n_2)$	$X_{(1)}^{(2)},\cdots,X_{(n_2)}^{(2)}$	$\bar{X}^{(2)}$
$\cdots$	$\cdots$	$\cdots$
$X_k(p\times n_k)$	$X_{(1)}^{(k)},\cdots,X_{(n_k)}^{(k)}$	$\bar{X}^{(k)}$

具体步骤：

• 作统计假设：$H_0:\mu_1=\cdots =\mu_k,H_1:至少存在i\neq j使得\mu_i\neq \mu_j$
• 计算离差阵: $$\begin{cases} 总离差阵：S_T=\sum_{i=1}^{k}\sum_{j=1}^{n_j}(X_{(j)}^{(i)}-\bar{X})(X_{(j)}^{(i)}-\bar{X})'\\ 组间离差阵：S_A=\sum_{i=1}^{k}n_i(\bar{X}^{(i)}-\bar{X})(\bar{X}^{(i)}-\bar{X})'\\ 组内离差阵：S_E=\sum_{i=1}^{k}\sum_{j=1}^{n_j}(X_{(j)}^{(i)}-\bar{X}^{(i)})(X_{(j)}^{(i)}-\bar{X}^{(i)})'\\ 离差阵分解公式：S_T=S_A+S_E\\ \end{cases}$$
• 计算统计量的具体值:
  $$\Lambda =\frac{\left |S_E\right |}{\left |S_E+S_A\right |}=\frac{\left |S_E\right |}{\left |S_T\right |}\sim \Lambda(p,n-k,k-1)$$
  $$-rln{\Lambda}\sim \chi^2(p(k-1)),r=(n-k)-\frac{1}{2}(p-(k-1)+1)$$
  其中$S_A\sim W_p(k-1,\sum),S_E\sim W_p(n-k,\sum)$
• 按规定的小概率标准$\alpha$，查$wilks$分布表并做出判断：
  -当${\Lambda}< \Lambda_{1-\alpha}(p,n-k.k-1)$，拒绝原假设，即认为有显著差异

*R实现：

#### k independent normal distribution ######### 
###################################################################
## H0: mu1=mu2=...=muk
## this is asymptotically a Chisq testing
##############  Input  ############################################
## data  = design matrix with a group index ind
############## Output  ############################################
## p.value     = p value
###################################################################
multi.mu.test=function(data, k)            

{
  ind=data$ind

  n=nrow(data)
  p=ncol(data)-1

  data=data[ ,1:p]
  T=(n-1)*var(data)

  A=0
  for (i in 1:k)                                
  {
    datai=data[ind==i, ]
    ni=nrow(datai)                                 
    A=A+(ni-1)*var(datai)
  }

  Lambda=det(A)/det(T)
  n1=n-k
  n2=k-1
  r=n1-(p-n2+1)/2
  Chi=(-1)*r*log(Lambda)

  p.value=1-pchisq(Chi, p*n2)
  return(p.value=p.value)
}