多元统计正态性检验作业
3.13 （1）对每个分量检验是否是一元正态分布 1.一维边缘分布的正态性检验 Q-Q图检验法
>data1=data.frame(x1=c(260,200,240,170,270,205,190,200,250,200,225,210,170,270,190,280,310,270,250,260),
x2=c(75,72,87,65,110,130,69,46,117,107,130,125,64,76,60,81,119,57,67,135),x3=c(40,34,45,39,39,34,27,45,21,28,36,26,31,33,34,20,25,31,31,39), x4=c(18,17,18,17,24,23,15,15,20,20,11,17,14,13,16,18,15,8,14,29)) >data2=data.frame(x1=c(310,310,190,225,170,210,280,210,280,200,200,280,190,295,270,280,240,280,370,280),x2=c(122,60,40,65,65,82,67,38,65,76,76,94,60,55,125,120,62,69,70,40),
x3=c(30,35,27,34,37,31,37,36,30,40,39,26,33,30,24,32,32,29,30,37), + x4=c(21,18,15,16,16,17,18,17,23,17,20,11,17,16,21,18,20,20,20,17)) >data3=data.frame(x1=c(320,260,360,295,270,380,240,260,260,295,240,310,330,345,250,260,225,345,360,250),x2=c(64,59,88,100,65,114,55,55,110,73,114,103,112,127,62,59,100,120,107,117),x3=c(39,37,28,36,32,36,42,34,29,33,38,32,21,24,22,21,34,36,25,36),x4=c(17,11,26,12,21,21,10,20,20,21,18,18,11,20,16,19,30,18,23,16)) > data=rbind(data1,data2,data3)
> qqnorm(data[,1]);qqline(data1[,1])
Normal Q-Q PlotSample Quantiles200250300350-2-10Theoretical Quantiles12
> qqnorm(data[,2]);qqline(data1[,2])
Normal Q-Q PlotSample Quantiles406080100120-2-10Theoretical Quantiles12
> qqnorm(data[,3]);qqline(data1[,3])
Normal Q-Q PlotSample Quantiles202530354045-2-10Theoretical Quantiles12
> qqnorm(data[,4]);qqline(data1[,4])
Normal Q-Q PlotSample Quantiles1015202530-2-10Theoretical Quantiles12
2.二元数据的正态性检验 1.等椭圆检验法
以检验（X1，X2）是否服从二元正态分布为例
> datax1x2=as.matrix(cbind(data[,1],data[,2])) > mean1=apply(datax1x2,2,mean);mean1 [1] 259.08333 84.11667 > s1=cov(datax1x2);s1 [,1] [,2] [1,] 2787.7048 433.6681 [2,] 433.6681 753.9692 > D2=c()
> for(i in 1:n){
+ D2[i]=(datax1x2[i,]-mean1)%*%solve(s1)%*%t(t(datax1x2[i,]-mean1)) + cat(D2[i])} > D2
[1] 0.1251861 1.2646933 0.1805204 2.8870075 0.8947723 5.3436674 1.7397987 2.4709765 1.7429712 2.7508022 4.2330751 4.2933030 2.9037836 [14] 0.1830727 1.9723611 0.2160584 1.9889309 1.2522066 0.3888125 3.7506293 2.2376907 2.4250881 3.3340246 0.6947016 2.8870075 0.9085821 [27] 0.7613934 3.0212101 0.8859113 1.2539095 1.2539095 0.2209607 1.9723611 2.2170959 2.2793965 1.7078104 0.6647098 0.6485288 5.8467138 [40] 3.4256245 2.6068141 0.9296899 3.8566119 0.6171641 0.6738596 5.4233607 1.1265728 1.2473480 0.9654576 0.8691904 1.7027351 1.1050343 [53] 2.2176911 3.9171163 0.6539759 0.9296899 1.0704952 3.3865331 3.7284302 1.7429712 > po=0.5
> p=dim(data)[2];p [1] 4
> d0=qchisq(p0,p);d0 [1] 3.356694
> a=sum(D2 > pi=a/n;pi [1] 0.8166667 #p0取0.5时，马氏距离小于d0的个数为49，占总样品比例约为0.82,拒绝来自二元正态分布的假设 > p0=0.25 > d0=qchisq(p0,p);d0 [1] 1.922558 > a=sum(D2 [1] 0.55 #p0取0.5时，马氏距离小于d0的个数为33，占总样品比例约为0.55,拒绝来自二元正态分布的假设 2.二元数据的?2图检验法 > pt=c() > for(t in 1:n){ + pt[t]=(t-0.5)/n + cat(pt[t])} > pt #pt [1] 0.008333333 0.025000000 0.041666667 0.058333333 0.075000000 0.091666667 0.108333333 0.125000000 0.141666667 0.158333333 0.175000000 [12] 0.191666667 0.208333333 0.225000000 0.241666667 0.258333333 0.275000000 0.291666667 0.308333333 0.325000000 0.341666667 0.358333333 [23] 0.375000000 0.391666667 0.408333333 0.425000000 0.441666667 0.458333333 0.475000000 0.491666667 0.508333333 0.525000000 0.541666667 [34] 0.558333333 0.575000000 0.591666667 0.608333333 0.625000000 0.641666667 0.658333333 0.675000000 0.691666667 0.708333333 0.725000000 [45] 0.741666667 0.758333333 0.775000000 0.791666667 0.808333333 0.825000000 0.841666667 0.858333333 0.875000000 0.891666667 0.908333333 [56] 0.925000000 0.941666667 0.958333333 0.975000000 0.991666667 > D2t=sort(D2) #D2(t) > xt2=c() > for(t in 1:n){ + xt2[t]=qchisq(pt[t],p) + cat(xt2[t]) + } > xt2 #?t2 [1] 0.2700151 0.4844186 0.6415772 0.7757695 0.8969359 1.0096230 1.1163677 1.2187621 1.3178880 1.4145247 1.5092595 1.6025523 [13] 1.6947743 1.7862337 1.8771930 1.9678806 2.0584996 2.1492342 2.2402545 2.3317204 2.4237845 2.5165951 2.6102978 2.7050379 [25] 2.8009620 2.8982198 2.9969656 3.0973602 3.1995725 3.3037815 3.4101784 3.5189686 3.6303748 3.7446397 3.8620297 3.9828389 [37] 4.1073944 4.2360619 4.3692534 4.5074361 4.6511434 4.8009895 4.9576873 5.1220712 5.2951282 5.4780385 5.6722300 5.8794549 [49] 6.1018972 6.3423292 6.6043460 6.8927308 7.2140471 7.5776562 7.9975859 8.4962822 9.1131220 9.9275079 11.1432868 13.6954281 > plot(D2t,pt) pt0.000.20.40.60.81.0123D2t456 （2）?2图检验对三组观测数据分别检验是否来自4元正态分布 对（1）组： > s1=cov(data1) > n1=dim(data1)[1] > mean1=apply(data1,2,mean) > data10=as.matrix(data1) > D2=c() > for(i in 1:n1){ + D2[i]=(data10[i,]-mean1)%*%solve(s1)%*%t(t(data10[i,]-mean1)) + cat(D2[i])} 2.3566150.87569193.3047952.8114523.7483283.4172392.5699034.3461183.5919072.10211511.080623.6419633.0200982.0029151.7946376.334355.2382985.1832451.2706137.3091 > D2t=sort(D2) #D2(t) > pt=c() > for(t in 1:n1){ + pt[t]=(t-0.5)/n1 + cat(pt[t])} #pt 0.0250.0750.1250.1750.2250.2750.3250.3750.4250.4750.5250.5750.6250.6750.7250.7750.8250.8750.9250.975 > plot(D2t,pt)