多元正态分布的协方差检验.docVIP

一、One sample covariance test cov.equal=function(x,Sigma,a=0.05) { ## x is the data set ## Sigma is the assumed covariance matrix ## a is the significance level set by default to 0.05 x=as.matrix(x) Sigma=as.matrix(Sigma) p=ncol(x) ## dimensionality of the data n=nrow(x) ## total sample size S=cov(x) ## sample covariance matrix## the next 2 lines construct the test statistic? mesa=solve(Sigma)%*%S test=sum(diag(mesa))-n*log(det(mesa))-n*p+n*p*log(n) df=0.5*p*(p+1) ## the degrees of freedom of the chi-square distribution pvalue=1-pchisq(test,df) ## p-value of the test statistic crit=qchisq(1-a,df) ## critical value of the chi-square distribution list(test=test,degres=df,p.value=pvalue,critical=crit) } 二、Multi-sample covariance matrices Log-likelihood ratio test cov.likel=function(x,ina,a=0.05) { ## x is the data set ## ina is a numeric vector indicating the groups of the data set ## a is the significance level, set to 0.05 by default x=as.matrix(x) p=ncol(x) ## dimension of the data set n=nrow(x) ## total sample size k=max(ina) ## number of groups nu=rep(0,k) ## the sample size of each group will be stored later pame=rep(0,k) ## the next 2 for functions separate the k groups and extract the ## covariance matrix of each group ## the way is not the best but it works nu=as.vector(table(ina))? mat=mat1=array(dim=c(p,p,k)) ## the next 3 lines create the pooled covariance matrix ## and calculate the covariance matrix of each group for (i in 1:k) {? mat[,,i]=((nu[i]-1)/nu[i])*cov(x[ina==i,])? mat1[,,i]=(nu[i]-1)*cov(x[ina==i,]) } Sp=apply(mat1,1:2,sum)/n for(iin1:k) pame[i]=det(solve(mat[,,i])%*%Sp) test=sum(nu*log(pame)) ## test statistic df=0.5*p*(p+1)*(k-1) ## degrees of freedom of the asymptotic chi-square pvalue=1-pchisq(test,df) ## p-value of the test statistic crit=qchisq(1-a,df) ## critical value of the chi-square distribution list(test=test,degrees=df,critical=crit,p.value=pvalue) Box’s M test cov.Mtest=function(x,ina,a=0.05) { ## x is the data set ## ina is a nume