《多元正态分布的协方差检验》由会员分享，可在线阅读，更多相关《多元正态分布的协方差检验（3页珍藏版）》请在金锄头文库上搜索。
1、一 、 One sample covariance testcov.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.05x=as.matrix(x)Sigma=as.matrix(Sigma)p=ncol(x) # dimensionality of the datan=nrow(x) # total sample sizeS=cov(x) # sample cov
2、ariance matrix# the next 2 lines construct the test statistic mesa=solve(Sigma)%*%Stest=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 distributionpvalue=1-pchisq(test,df) # p-value of the test statisticcrit=qchisq(1-a,df) # critical value of t
3、he chi-square distributionlist(test=test,degres=df,p.value=pvalue,critical=crit) 二 、 Multi-sample covariance matricesLog-likelihood ratio testcov.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
4、by defaultx=as.matrix(x)p=ncol(x) # dimension of the data setn=nrow(x) # total sample sizek=max(ina) # number of groupsnu=rep(0,k) # the sample size of each group will be stored laterpame=rep(0,k)# the next 2 for functions separate the k groups and extract the# covariance matrix of each group# the w
5、ay is not the best but it worksnu=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 groupfor (i in 1:k) mat,i=(nui-1)/nui)*cov(xina=i,) mat1,i=(nui-1)*cov(xina=i,) Sp=apply(mat1,1:2,sum)/nfor(iin1:k)pam
6、ei=det(solve(mat,i)%*%Sp)test=sum(nu*log(pame) # test statisticdf=0.5*p*(p+1)*(k-1) # degrees of freedom of the asymptotic chi-squarepvalue=1-pchisq(test,df) # p-value of the test statisticcrit=qchisq(1-a,df) # critical value of the chi-square distributionlist(test=test,degrees=df,critical=crit,p.va
7、lue=pvalue) Boxs M testcov.Mtest=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 defaultx=as.matrix(x)p=ncol(x) # dimension of the data setn=nrow(x) # total sample sizek=max(ina) # number of groups
8、nu=rep(0,k) # the sample size of each group will be stored here laterpame=rep(0,k) # the determinant of each covariance will be stored here# the next for function calculates the covariance matrix of each groupnu=as.vector(table(ina) mat=mat1=array(dim=c(p,p,k)for (i in 1:k) mat,i=cov(xina=i,)pamei=d
9、et(mat,i) # the detemirnant of each covariance matrix mat1,i=(nui-1)*cov(xina=i,) # the next 2 lines calculate the pooled covariance matrixSp=apply(mat1,1:2,sum)Sp=Sp/(n-k)for (i in 1:k)pamela=det(Sp) # determinant of the pooled covariance matrixtest1=sum(nu-1)*log(pamela/pame)gama1=(2*(p2)+3*p-1)/(
10、6*(p+1)*(k-1)gama2=(sum(1/(nu-1)-1/(n-k)gama=1-gama1*gama2test=gama*test1 # this is the M (test statistic)df=0.5*p*(p+1)*(k-1) # degrees of freedom of the chi-square distributionpvalue=1-pchisq(test,df) # p-value of the test statisticcrit=qchisq(1-a,df) # critical value of the chi-square distributionlist(M.test=test,degrees=df,critical=crit,p.value=pvalue)
