本文關(guān)鍵詞：多元統(tǒng)計(jì)分析課程，由筆耕文化傳播整理發(fā)布。

引言

本學(xué)期也開了一門多元統(tǒng)計(jì)分析課程，也趁機(jī)想把課后上機(jī)題實(shí)現(xiàn)一遍，以增強(qiáng)理解。

教材使用的是約翰遜的《多元統(tǒng)計(jì)分析》第六版，中英文版教材、數(shù)據(jù)集、講義見
還參考了王斌會(huì)老師的《多元統(tǒng)計(jì)分析及R語(yǔ)言建�！�

本文內(nèi)容主要為第4章多元正態(tài)分布的上機(jī)題，圖略。
[rmd文檔見](）
可以直接用Rstudio打開（之前先安裝knitr包）

4.28 data_4.28<-read.table("E:\\研究生\\應(yīng)用多元統(tǒng)計(jì)\\JohnsonWichern Data sets\\T1-5.DAT") #正態(tài)Q-Q圖 qqnorm(data_4.28$V2) #正態(tài)性檢驗(yàn) #原始數(shù)據(jù)排序 new_data<-sort(data_4.28$V2) length(new_data) #對(duì)應(yīng)概率值 prob<-((i-0.5)/n) } all_pro<-sapply(1:42,prob)#所有概率值 #對(duì)應(yīng)的標(biāo)準(zhǔn)正態(tài)分位數(shù) all_q<-qnorm(all_pro) #Q-Q圖的相關(guān)系數(shù) rq<-cor(new_data,all_q) #由于Q-Q圖的相關(guān)系數(shù)rq為0.9693258，小于表4-2中n=40對(duì)應(yīng)的臨界點(diǎn)，，所以拒絕正態(tài)性假設(shè)。 4.29 #(a) #計(jì)算樣本協(xié)方差矩陣 s<-cov(data_4.28[,5:6]) #s的逆 s_solve<-solve(s) x_bar<-apply(data_4.28[,5:6],MARGIN=2,mean)#兩列平均數(shù) x_bar<-matrix(as.vector(x_bar),42,2,by=2) two_col<-t(data_4.28[,5:6]-x_bar)#兩列x-x_bar #計(jì)算所用統(tǒng)計(jì)距離dis dis<-c() for(i in 1:length(two_col[1,])){ dis[i]<-t(two_col[,i])%*%s_solve%*%two_col[,i] } chisq_num<-qchisq(0.5,2) #所占比例 pro<-length(which(dis<chisq_num))/length(dis) sort_data<-sort(dis) #概率密度為4.28中的all_pro #對(duì)應(yīng)的自由度為2的卡方分位數(shù) all_chiisq<-sapply(all_pro,qchisq,df=2)#所有概率值 #畫出卡方圖 也就是(all_chiisq,sort_data)對(duì)應(yīng)的散點(diǎn)圖 library(ggplot2) qplot(all_chiisq, sort_data, geom='point') 4.30 #讀入數(shù)據(jù) 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) #構(gòu)建冪變化函數(shù) ##冪類變化函數(shù)(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) } #生成多個(gè)λ，求使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))] #進(jìn)行數(shù)據(jù)冪變化 new_data<-sapply(data_4.30_x1,box_cox,λ=max_λ) #變化后的Q-Q圖 qqnorm(new_data) λ<-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))] #進(jìn)行數(shù)據(jù)冪變化 new_data<-sapply(data_4.30_x2,box_cox,λ=max_λ) #變化后的Q-Q圖 qqnorm(new_data) 4.39 data_4.] norm_test<-function(data){ #原始數(shù)據(jù)排序 new_data<-sort(data) len_data<-length(new_data) prob<-function(i,n){#構(gòu)建一個(gè)概率值的函數(shù) return((i-0.5)/n) } #對(duì)應(yīng)概率值 all_pro<-sapply(all_q<-qnorm(all_pro) #Q-Q圖的相關(guān)系數(shù) return(cor(new_data,all_q)) } ##對(duì)于獨(dú)立性 #Q-Q圖 qqnorm(data_4.39$V1)#大部分在一條直線上 norm_test(data_4.39$V1) qqnorm(data_4.39$V2)#大部分在一條直線上 norm_test(data_4.39$V2) #在顯著性水平為0.05的情況下，當(dāng)n=150時(shí)，0.989小于表4.2中的0.9913拒絕正態(tài)性假定 ##對(duì)于仁愛心 qqnorm(data_4.39$V3)#大部分在一條直線上 norm_test(data_4.39$V3) #在顯著性水平為0.05的情況下，當(dāng)n=150時(shí)，0.993大于表4.2中的0.9913不拒絕正態(tài)性假定 #對(duì)于順從性 qqnorm(data_4.39$V4)#大部分在一條直線上 norm_test(data_4.39$V4) #在顯著性水平為0.05的情況下，當(dāng)n=150時(shí)，0.993大于表4.2中的0.9913 不拒絕正態(tài)性假定 #對(duì)于領(lǐng)導(dǎo)能力 qqnorm(data_4.39$V5)#大部分在一條直線上 norm_test(data_4.39$V5) chis_chart<-function(x){ #計(jì)算樣本協(xié)方差矩陣 s<-cov(x) #s的逆 s_solve<-solve(s) x_bar<-apply(x,MARGIN=2,mean)#兩列平均數(shù) two_col<-t(x-x_bar)#兩列x-x_bar #計(jì)算所用統(tǒng)計(jì)距離dis dis<-c() (two_col[1,])){ dis[i]<-t(two_col[,i])%*%s_solve%*%two_col[,i] } #對(duì)廣義平方距離dis進(jìn)行排序 sort_data<-sort(dis) #prob在題4.28中構(gòu)造 all_pro<-sapply(all_chiisq<-sapply(all_pro,qchisq,df=5)#所有概率值 #畫出卡方圖 也就是(all_chiisq,sort_data)對(duì)應(yīng)的散點(diǎn)圖 library(ggplot2) qplot(all_chiisq, sort_data, geom='point') } chis_chart(data_4.39) λ<-seq(-1,2,0.1) all_l<-c() (λ)){ all_l[n]<-l_value(data_4.39$V1,lamda=λ[n]) } #取使變化后的l_value最大的λ值 max_λ<-λ[which(all_l==max(all_l))] #進(jìn)行數(shù)據(jù)冪變化 new_data<-sapply(data_4.39$V1,box_cox,λ=max_λ) #變化后的Q-Q圖 qqnorm(new_data) ##對(duì)于支撐力 all_l<-c() (λ)){ all_l[n]<-l_value(data_4.39$V2,lamda=λ[n]) } #取使變化后的l_value最大的λ值 max_λ<-λ[which(all_l==max(all_l))] #進(jìn)行數(shù)據(jù)冪變化 new_data<-sapply(data_4.39$V2,box_cox,λ=max_λ) #變化后的Q-Q圖 qqnorm(new_data) ##對(duì)于領(lǐng)導(dǎo)力 all_l<-c() (λ)){ all_l[n]<-l_value(data_4.39$V5,lamda=λ[n]) } #取使變化后的l_value最大的λ值 max_λ<-λ[which(all_l==max(all_l))] #進(jìn)行數(shù)據(jù)冪變化 new_data<-sapply(data_4.39$V5,box_cox,λ=max_λ) #變化后的Q-Q圖 qqnorm(new_data) 4.40 data_4.40<-read.table("E:\\研究生\\應(yīng)用多元統(tǒng)計(jì)\\JohnsonWichern Data sets\\T1-11.DAT") library(ggplot2) #散點(diǎn)圖檢查 qplot(data_4.40$V1, data_4.40$V2, geom='point') #從散點(diǎn)圖可以看出在x軸和y軸分別有一個(gè)離群值 #標(biāo)準(zhǔn)化值來檢查 cen_data<-scale(data_4.40) #每一列的最大離群值為 apply(abs(cen_data),2,max) #與取標(biāo)準(zhǔn)化數(shù)據(jù)比較，第一列第13行，第二列第7行與其他數(shù)據(jù)存在較大偏離 #(b)(c)略4.40略

  本文關(guān)鍵詞：多元統(tǒng)計(jì)分析課程，由筆耕文化傳播整理發(fā)布。

本文編號(hào)：195424