多元正态分布的协方差检验、独立性检验及R实现多元正态分布的协方差检验、

1.单个 pp 元正态总体协方差阵的检验具体步骤：

作统计假设(1)： H0:∑=Ip,H1:∑≠IpH 0 : ∑ = I p , H 1 : ∑ ≠ I p
统计量： λ=supθ∈Θ0L(μ,Ip)supθ∈ΘL(μ,∑)=(en)np/2|S|n/2exp[?12tr(S)]λ = s u p θ ∈ Θ 0 L ( μ , I p ) s u p θ ∈ Θ L ( μ , ∑ ) = ( e n ) n p / 2 | S | n / 2 e x p [ ? 1 2 t r ( S ) ]Λ=?2lnλ～χ2(p(p+1)2)Λ = ? 2 l n λ ～ χ 2 ( p ( p + 1 ) 2 ) 其中 supθ∈Θ0L(μ,Ip)=L(Xˉ,Ip)=(2π)?np/2exp[?12tr(S)]s u p θ ∈ Θ 0 L ( μ , I p ) = L ( X ˉ , I p ) = ( 2 π ) ? n p / 2 e x p [ ? 1 2 t r ( S ) ]supθ∈ΘL(μ,∑)=L(Xˉ,Sn)=(2π)?np/2∣∣∣1nS∣∣∣?n/2exp[?12pn]s u p θ ∈ Θ L ( μ , ∑ ) = L ( X ˉ , S n ) = ( 2 π ) ? n p / 2 | 1 n S | ? n / 2 e x p [ ? 1 2 p n ]
作统计假设(2)： H0:∑=∑0,H1:∑≠∑0H 0 : ∑ = ∑ 0 , H 1 : ∑ ≠ ∑ 0
统计量：做变换转换为上面的假设： Yi=CXi,CY i = C X i , C 为非奇异矩阵
λ=(en)np/2|SY|n/2exp[?12tr(SY)]λ = ( e n ) n p / 2 | S Y | n / 2 e x p [ ? 1 2 t r ( S Y ) ] 其中 SY=CS′CS Y = C S ′ C
作统计假设(3)：（球性检验） H0:∑=σ2∑0H 0 : ∑ = σ 2 ∑ 0 ( σ2未知σ 2 未知 )
统计量：
σ2σ 2 的极大似然估计为 1nptr(∑?10S)1 n p t r ( ∑ 0 ? 1 S )
λ=∣∣∑?10S∣∣n/2[tr(∑?10S)/p]np/2λ = | ∑ 0 ? 1 S | n / 2 [ t r ( ∑ 0 ? 1 S ) / p ] n p / 2W=(λ)2/n=pp∣∣∑?10S∣∣[tr(∑?10S)]pW = ( λ ) 2 / n = p p | ∑ 0 ? 1 S | [ t r ( ∑ 0 ? 1 S ) ] p?((n?1)?2p2+p+26p)lnW～χ2(p(p+1)2?1)? ( ( n ? 1 ) ? 2 p 2 + p + 2 6 p ) l n W ～ χ 2 ( p ( p + 1 ) 2 ? 1 )

*R实现：

#### variance testing ######### var.test=function(data, Sigma0) ############################################################### ## H0: Sigma=Sigma0 ## this is aymptotically a Chisq testing ##############Input######################################## ## data= https://www.it610.com/article/design matrix with the ith sample in the ith line ## Sigma0= Simga0 for null hypothesis ############## Output######################################## ## p.value= p value ############################################################### { n=nrow(data) p=ncol(data)A=(n-1)*var(data) S=A%*%solve(Sigma0)lambda=exp(sum(diag((-1)*S/2)))*(det(S))^(n/2)*(exp(1)/n)^(n*p/2) T5=-2*log(lambda)p.value=1-pchisq(T5, p*(p+1)/2) return(p.value) }

2.两个 pp元正态总体协方差阵相等的检验具体步骤：

作统计假设： H0:∑1=∑2,H1:∑1≠∑2H 0 : ∑ 1 = ∑ 2 , H 1 : ∑ 1 ≠ ∑ 2
统计量：（式子太长有时间补上）

*R实现：

#### k independent normal distribution######### multi.var.test=function(data, k) ################################################################### ## H0: Sigma1=Sigma2=...=Sigmak ## this is asymptotically a Chisq testing ##############Input############################################ ## data= https://www.it610.com/article/design matrix with a group index ind ############## Output############################################ ## p.value= p value ################################################################### { ind=data$indn=nrow(data) p=ncol(data)-1 data=data[ ,1:p]A=0 for (i in 1:k) { datai=data[ind==i, ] ni=nrow(datai) A=A+(ni-1)*var(datai) }det.A=0 for (i in 1:k) { datai=data[ind==i, ] ni=nrow(datai) det.A=det.A+(ni-1)*log(det(var(datai))) }M=(n-k)*log(det(A/(n-k)))-det.A d=(2*p^2+3*p-1)*(k+1)/(6*(p+1)*(n-k)) f=p*(p+1)*(k-1)/2T6=(1-d)*Mp.value=1-pchisq(T6, f) return(p.value=p.value) }

3.独立性检验（二个总体） X～Np(μ,∑),X={X1X2},X1:q,X2:p?qX ～ N p ( μ , ∑ ) , X = { X 1 X 2 } , X 1 : q , X 2 : p ? q
具体步骤：

作统计假设： H0:∑12=0H 0 : ∑ 12 = 0
统计量：记∑0={∑1100∑22}记 ∑ 0 = { ∑ 11 0 0 ∑ 22 }S=∑i=1n(Xi?Xˉ)(Xi?Xˉ)′={S11S21S12S22}S = ∑ i = 1 n ( X i ? X ˉ ) ( X i ? X ˉ ) ′ = { S 11 S 12 S 21 S 22 }λ2/n=|S||S11||S22|～Λ(p?q)λ 2 / n = | S | | S 11 | | S 22 | ～ Λ ( p ? q )
对K个总体： V=|S|∏ki=1|Sii|=∏i=1k?1vi,vi～Λ(mi,n?∑j=1imj?1,∑j=1imj)V = | S | ∏ i = 1 k | S i i | = ∏ i = 1 k ? 1 v i , v i ～ Λ ( m i , n ? ∑ j = 1 i m j ? 1 , ∑ j = 1 i m j )?blnV～χ2(f)? b l n V ～ χ 2 ( f ) 其中 b=n?32?p3?∑k1p3i3(p2?∑k1p2i)b = n ? 3 2 ? p 3 ? ∑ 1 k p i 3 3 ( p 2 ? ∑ 1 k p i 2 )f=12[p(p+1)?∑1kpi(pi+1)]f = 1 2 [ p ( p + 1 ) ? ∑ 1 k p i ( p i + 1 ) ]

【多元正态分布的协方差检验、独立性检验及R实现】*R实现：

######### testing for independent######### multi.ind.test=function(data, k) ################################################################### ## H0: sigma(ij)=0,i!=j ## this is asymptotically a Chisq testing ##############Input############################################ ## data= https://www.it610.com/article/design matrix with a group index ind ############## Output############################################ ## p.value= p value ################################################################### { ind=data$ind #不同总体分组#n=nrow(data)-1 #因为最后一行为ind# p=ncol(data) data=data[1:n , ] X.bar=apply(data, 2, mean) S=(n-1)*var(data)det.B=1 P1=0 P2=0 P3=0 for (i in 1:k) { datai=data[ ,ind==i] p_i=ncol(datai) P1=P1+p_i**3 P2=P2+p_i**2 P3=P3+p_i*(p_i+1) Si=(n-1)*var(datai) det.B=det.Si*det.B }V=det.S/det.B b=n-3/2-(p^3-P1)/(3*p^2-3*P2) f=0.5*(p(p+1)-P3)T6=(-b)*log(V)p.value=1-pchisq(T6, f) return(p.value=p.value) }