多元统计分析上机题之R语言实现（多元正态分布）

引言

本学期也开了一门多元统计分析课程，也趁机想把课后上机题实现一遍，以增强理解。

4.28

data_4.28<-read.table("E:\\研究生\\应用多元统计\\JohnsonWichern Data sets\\T1-5.DAT") #正态Q-Q图 norm(data_4.28$V2) #正态性检验 #原始数据排序 new_data<-sort(data_4.28$V2) length(new_data) #对应概率值 prob<-function(i,n=42){ 
   #构建一个概率值的函数 return((i-0.5)/n) } all_pro<-sapply(1:42,prob)#所有概率值 #对应的标准正态分位数 all_q<-qnorm(all_pro) #Q-Q图的相关系数 rq<-cor(new_data,all_q) #由于Q-Q图的相关系数rq为0.，小于表4-2中n=40对应的临界点，所以拒绝正态性假设。

4.29

#(a) #计算样本协方差矩阵 s<-cov(data_4.28[,5:6]) #s的逆 s_solve<-solve(s) x_bar<-apply(data_4.28[,5:6],MARGIN=2,mean)#两列平均数 x_bar<-matrix(as.vector(x_bar),42,2,by=2) two_col<-t(data_4.28[,5:6]-x_bar)#两列x-x_bar #计算所用统计距离dis dis<-c() for(i in 1:length(two_col[1,])){ dis[i]<-t(two_col[,i])%*%s_solve%*%two_col[,i] }  #(b) #自由度为2概率密度为0.5的卡方分布临界值 chisq_num<-qchisq(0.5,2) #所占比例 pro<-length(which(dis 
   
    
    
    
    
    
   #(c) 
   #对广义平方距离dis进行排序 sort_data<-sort(dis) 
   #概率密度为4.28中的all_pro 
   #对应的自由度为2的卡方分位数 all_chiisq<-sapply(all_pro,qchisq,df= 
   2) 
   #所有概率值 
   #画出卡方图 也就是(all_chiisq,sort_data)对应的散点图 library(ggplot2) qplot(all_chiisq, sort_data, geom= 
   'point') 
  

4.30

#读入数据 data_4.30_x1<-c(1:9,11) data_4.30_x2<-c(18.95,19.00,17.95,15.54,14.00,12.95,8.94,7.49,6.00,3.99) #构建幂变化函数 幂类变化函数(Box-Cox) box_cox<-function (x,λ){ if (λ==0) { return(log(x)) }else{ return((x^λ-1)/λ) } } l_value<-function(X,lamda){ x_new<-sapply(X,box_cox,λ=lamda) x_bar<-mean(x_new) l_val<-log(mean((x_new-x_bar)^2))*(-length(x_new)/2)+(lamda-1)*sum(log(X)) return(l_val) } #生成多个λ，求使l_value最大的λ_hat值 λ<-seq(-1,2,0.1) all_l<-c() for(n in 1:length(λ)){ all_l[n]<-l_value(data_4.30_x1,lamda=λ[n]) } #取使变化后的l_value最大的λ值 max_λ<-λ[which(all_l==max(all_l))] #进行数据幂变化 new_data<-sapply(data_4.30_x1,box_cox,λ=max_λ) #变化后的Q-Q图 norm(new_data)  #(b) #基本同(a)题 λ<-seq(-1,2,0.1) all_l<-c() for(n in 1:length(λ)){ all_l[n]<-l_value(data_4.30_x2,lamda=λ[n]) } #取使变化后的l_value最大的λ值 max_λ<-λ[which(all_l==max(all_l))] #进行数据幂变化 new_data<-sapply(data_4.30_x2,box_cox,λ=max_λ) #变化后的Q-Q图 norm(new_data)  #(c)略 #题4.31-4.38均按照4.28-4.30的解题思路进行即 #考虑边缘正态性：先做Q-Q图做个粗略的了解 然后计算Q-Q图的相关系数 并与书中表4.2进行比较 得出是否拒绝正态性的假设 #考虑二维正态性 采用4.29的方法 做卡方图 #变换可以采用平方根变换 对数变换 z变换 ，见书本p147页，还可以使用4.30中的幂变换，然后将变换后的数据画Q-Q图进行判断。

4.39

data_4.39<-read.table("E:\\研究生\\应用多元统计\\JohnsonWichern Data sets\\T4-6.DAT")[,1:5] #(a) #正态性检验 #计算Q-Q图的相关系数 将题4.28的代码进行封装 norm_test<-function(data){ #原始数据排序 new_data<-sort(data) len_data<-length(new_data) prob<-function(i,n){ 
  #构建一个概率值的函数 return((i-0.5)/n) } #对应概率值 all_pro<-sapply(1:len_data,prob,n=len_data)#所有概率值 #对应的标准正态分位数 all_q<-qnorm(all_pro) #Q-Q图的相关系数 return(cor(new_data,all_q)) } 对于独立性 #Q-Q图 norm(data_4.39$V1)#大部分在一条直线上 norm_test(data_4.39$V1) #在显著性水平为0.05的情况下，当n=150时，0.988小于于表4.2中的0.9913拒绝正态性假定。 #也可以采用shapiro-wilk检验 #使用在mvnormtest包里mshapiro.test，具体可以使用?mshapiro.test查看使用方法 对于支撑力 norm(data_4.39$V2)#大部分在一条直线上 norm_test(data_4.39$V2) #在显著性水平为0.05的情况下，当n=150时，0.989小于表4.2中的0.9913拒绝正态性假定 对于仁爱心 norm(data_4.39$V3)#大部分在一条直线上 norm_test(data_4.39$V3) #在显著性水平为0.05的情况下，当n=150时，0.993大于表4.2中的0.9913不拒绝正态性假定 #对于顺从性 norm(data_4.39$V4)#大部分在一条直线上 norm_test(data_4.39$V4) #在显著性水平为0.05的情况下，当n=150时，0.993大于表4.2中的0.9913 不拒绝正态性假定 #对于领导能力 norm(data_4.39$V5)#大部分在一条直线上 norm_test(data_4.39$V5) #在显著性水平为0.05的情况下，当n=150时，0.981小于表4.2中的0.9913 拒绝正态性假定  #(b) 使用卡方图进行判定 #构造画卡方图的函数 方法同题4.29 chis_chart<-function(x){ #计算样本协方差矩阵 s<-cov(x) #s的逆 s_solve<-solve(s) x_bar<-apply(x,MARGIN=2,mean)#两列平均数 two_col<-t(x-x_bar)#两列x-x_bar #计算所用统计距离dis dis<-c() for(i in 1:length(two_col[1,])){ dis[i]<-t(two_col[,i])%*%s_solve%*%two_col[,i] } #对广义平方距离dis进行排序 sort_data<-sort(dis) #prob在题4.28中构造 all_pro<-sapply(1:length(x[,1]),prob,n=130)#所有概率值 #对应的自由度为5的卡方分位数 all_chiisq<-sapply(all_pro,qchisq,df=5)#所有概率值 #画出卡方图 也就是(all_chiisq,sort_data)对应的散点图 library(ggplot2) qplot(all_chiisq, sort_data, geom='point') } chis_chart(data_4.39) #很明显，卡方图上点不是接近于一条直线，偏一条曲线，所以多元正态性不满足，可知，边缘正态性不满足的情况下，多元正态性也很少满足 #  #(c) #在(a)中，独立性、支撑力、领导力的分布不符合正态性 幂变化函数构造见题4.30 对于独立性 #生成多个λ，求使l_value最大的λ_hat值 λ<-seq(-1,2,0.1) all_l<-c() for(n in 1:length(λ)){ all_l[n]<-l_value(data_4.39$V1,lamda=λ[n]) } #取使变化后的l_value最大的λ值 max_λ<-λ[which(all_l==max(all_l))] #进行数据幂变化 new_data<-sapply(data_4.39$V1,box_cox,λ=max_λ) #变化后的Q-Q图 norm(new_data) 对于支撑力 all_l<-c() for(n in 1:length(λ)){ all_l[n]<-l_value(data_4.39$V2,lamda=λ[n]) } #取使变化后的l_value最大的λ值 max_λ<-λ[which(all_l==max(all_l))] #进行数据幂变化 new_data<-sapply(data_4.39$V2,box_cox,λ=max_λ) #变化后的Q-Q图 norm(new_data) 对于领导力 all_l<-c() for(n in 1:length(λ)){ all_l[n]<-l_value(data_4.39$V5,lamda=λ[n]) } #取使变化后的l_value最大的λ值 max_λ<-λ[which(all_l==max(all_l))] #进行数据幂变化 new_data<-sapply(data_4.39$V5,box_cox,λ=max_λ) #变化后的Q-Q图 norm(new_data) 

4.40

data_4.40<-read.table("E:\\研究生\\应用多元统计\\JohnsonWichern Data sets\\T1-11.DAT") library(ggplot2) #散点图检查 qplot(data_4.40$V1, data_4.40$V2, geom='point') #从散点图可以看出在x轴和y轴分别有一个离群值 #标准化值来检查 cen_data<-scale(data_4.40) #每一列的最大离群值为 apply(abs(cen_data),2,max) #与取标准化数据比较，第一列第13行，第二列第7行与其他数据存在较大偏离 #(b)(c)略4.40略