%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%应用多元统计分析 2010.04.02
%%%%正态性检验
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

> stiff<-read.csv("Stiff.csv")
> head(stiff)
    X1   X2   X3   X4  Dsq
1 1889 1651 1561 1778 0.60
2 2403 2048 2087 2197 5.48
3 2119 1700 1815 2222 7.62
4 1645 1627 1110 1533 5.21
5 1976 1916 1614 1883 1.40
6 1712 1712 1439 1546 2.22
> dim(stiff)
[1] 30  5

1. Pearson's chi-squared test

> A<-table(cut(stiff[,1],br=c(0,1500,2000,2500,3000)))
> p<-pnorm(c(1500,2000,2500,3000),mean(stiff[,1]),sd(stiff[,1]))
> p<-c(p[1],p[2]-p[1],p[3]-p[2],1-p[3])
> chisq.test(A,p=p)

        Chi-squared test for given probabilities

data:  A 
X-squared = 1.1368, df = 3, p-value = 0.7682

Warning message:
In chisq.test(A, p = p) : Chi-squared近似算法有可能不准

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2.K-S test
> ks.test(stiff[,1],"pnorm",mean(stiff[,1]),sd(stiff[,1]))

        One-sample Kolmogorov-Smirnov test

data:  stiff[, 1] 
D = 0.1149, p-value = 0.7819
alternative hypothesis: two-sided 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3. Shapiro-Wilk test
> shapiro.test(stiff[,1])

        Shapiro-Wilk normality test

data:  stiff[, 1] 
W = 0.9307, p-value = 0.05118

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4.Normal Q-Q 
> qqnorm(stiff[,1])
> qqline(stiff[,1])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%    多元正态性检验
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1.Chi-squared Q-Q plot

> distsq  #计算马氏距离
function(x)
{
        dif <- as.matrix(sweep(x, 2, 
                apply(x, 2, mean), "-"))
        diag(dif %*% solve(var(x)) %*% t(
                dif))
}

> qqchisq   #画Q-Q图
function(x, DF)    #x为马氏距离
{
        xv <- as.vector(x)
        plot(qchisq(ppoints(length(xv)), 
                DF), sort(xv), main="Chi-squared Q-Q Test",ylab = 
                "ordered data", xlab = 
                paste("Quantiles of Chi_Sq\n (",
                DF, "d.f.)"))
        abline(0, 1)
        invisible()
}

> dsq<-distsq(stiff[,1:4])
> qqchisq(dsq,4)

2.等概率椭圆检验法
> elliptest
function(x,alpha=0.5)
{
df<-dim(x)[2]
n<-dim(x)[1]
num<-0
 dif <- as.matrix(sweep(x, 2, 
                apply(x, 2, mean), "-"))
 dsq<- diag(dif %*% solve(var(x)) %*% t(
                dif))
for (i in 1:n) {
if (dsq[i]<=qchisq(alpha,df)) num<-num+1
}
ratio<-num/n
print(ratio) #输出样本落在椭圆内的比例
}

> elliptest(stiff[,1:4])
[1] 0.5666667
> elliptest(stiff[,1:4],alpha=0.7)
[1] 0.7
> elliptest(stiff[,1:4],alpha=0.8)
[1] 0.8333333
> elliptest(stiff[,1:4],alpha=0.6)
[1] 0.6666667