非參數統計分析

實驗内容及要求

1. 利用觀測資料計算總體分位數、對稱中心和位置差的點估計，區間估計；
2. 利用R軟體自帶程式或自程式設計式完成中位數的符号檢驗，兩總體比較的Wilcoxon秩和檢驗和K-S檢驗，獨立性與随機性的卡方檢驗和Fisher列聯表檢驗，相關性秩檢驗與協同性檢驗以及多總體比較的秩和檢驗和卡方檢驗；
3. 制作資料經驗分布函數、機率密度圖像，使用分布拟合方法解決總體類型的檢驗問題；
4. 通過最小二乘與權函數結合使用的方法解決非線性回歸問題；

本文中國運動員獲世界冠軍實驗項目資料來源于國家資料網站：

http://data.stats.gov.cn/easyquery.htm?cn=C01

實驗項目一

項目一收集1999年至2017年中國運動員獲世界冠軍的項數（機關：項）和人數（機關：人）。

年份 獲世界冠軍項數 獲世界冠軍人數 男子世界冠軍人數 女子世界冠軍人數 1999年 91 129 57 72 2000年 92 109 49 60 2001年 79 138 61 77 2002年 99 123 42 81 2003年 17 94 31 63 2004年 27 175 77 98 2005年 22 159 70 89 2006年 24 169 82 87 2007年 22 217 94 123 2008年 24 151 68 83 2009年 30 223 89 2010年 22 180 99 2011年 24 198 96 2012年 24 140 2013年 22 164 2014年 22 206 2015年 25 214 2016年 23 154 2017年 24 248

(1)計算其中運動員獲世界冠軍人數(人)的總體分位數如下：

$人數排序

94 109 123 129 138 140 151 154 159 164 169 175 180 198 206 214   217 223

248

$五數總括

94 139 164 202 248

$分位數

  0%  25%  50%  75% 100%

  94  139  164  202  248

(2) 計算其中運動員獲世界冠軍人數(人)的對稱中心的點估計、區間估計：

點估計：

$均值

[1] 37.52632

$截尾均值

[1] 167.9474

區間估計：求出置信度為0.9的置信區間

$資料位置

[1] "( 5 , 15 )"

$區間估計

[1] "[ 22 , 30 ]"

(3)設運動員獲世界冠軍項數(項)和運動員獲世界冠軍人數(人)這兩個簡單樣本分别取自總體Y和X，假定，試用樣本均值差作的估計。

樣本均值易受異常值影響，但是中位數不會，是以對于本樣本同時采用中位數來作為位置差的點估計。

$樣本均值之差

[1] 130.4211

$樣本中位數之差

[1] 140

位置差的HL區間估計： 

$中位數

[1] 131

$上下界位置

[1] "114" "248"

$`95%置信區間`

[1] "[ 107 , 152 ]

實驗項目二

(1) 中位數的符号檢驗

檢驗1999年至2017年期間，運動員獲世界冠軍人數(人)的中位數是否為155？()

中為負數的個數，則的拒絕域為

根據公式

由觀測值得到Y的值，若則接受，否則拒絕。

運動員獲世界冠軍人數(人)
P值	0.6476059

由于p值大于0.05，檢驗的結論是接受原假設，認為運動員獲世界冠軍人數(人)的中位數應該為155。

(2)兩總體比較的Wilcoxon秩和檢驗

運動員獲世界冠軍男女人數是否存在顯著差異？()

拒絕域為：

通過軟體計算出：

$r1

[1] 88

$r2

[1] 152

$是否落入拒絕域

[1] FALSE

則得到結論：運動員獲世界冠軍男女人數存在顯著差異。

(3) 兩總體比較的K-S檢驗

設男子世界冠軍人數來自分布為F(x)的總體女子世界冠軍人數來自分布為Gx)的總體，檢驗這兩個分布是否相同，即原假設為：

計算結果如下：

Two-sample Kolmogorov-Smirnov test

data:  x and y

D = 0.33846, p-value = 0.5366

alternative hypothesis: two-sided

由于p值=0.5366>0.05，故接受原假設，即認為F(x)和G(x)這兩個分布函數相同。

(4)卡方獨立性檢驗

本案例資料來源于百度文庫。

在遇到車禍的情況下，乘客系安全帶與沒系安全帶時受到的沖擊力的資料如下：

受傷情況	無	輕微	較重	嚴重	合計
安全帶系	12813	647	359	42	13861
安全帶沒系	65963	4000	2642	303	72908
合計	78776	4647	3001	345	86769

各因子之間是否是獨立的？

計算結果如下：

Pearson's Chi-squared test

data:  rbind(yesbelt, nobelt)

X-squared = 59.224, df = 3, p-value = 8.61e-13

其中p值=8.61e-13<0.05，拒絕原假設，即認為兩行因子不是互相獨立的。

(5) Fisher列聯表檢驗

本例仍然使用乘客系安全帶與沒系安全帶時受到的沖擊力的資料進行分析，首先将上例資料改為2*2的列聯表：

受傷情況	無	受傷	合計
安全帶系	12813	1048	13861
安全帶沒系	65963	6945	72908
合計	78776	7993	86769

計算結果如下：

Fisher's Exact Test for Count Data

data:  x

p-value = 6.971e-14

alternative hypothesis: true odds ratio is not equal to 1

95 percent confidence interval:

 1.202757 1.378612

sample estimates:

odds ratio

  1.287236

由于P值=6.971e-14 < 0.05，是以拒絕原假設，即認為受傷情況與是否系安全帶有顯著關系。由于賠率比大于1，是以還是正相關。

(6)Spearman秩相關檢驗方法

選取1999年到2008年運動員獲世界冠軍男女人數資料，檢驗二者之間是否存在相關關系？

計算結果如下：

Spearman's rank correlation rho

data:  x and y

S = 22, p-value = 0.002681

alternative hypothesis: true rho is not equal to 0

sample estimates:

      rho

0.8666667

結果顯示P值=0.002681 < 0.05，是以拒絕原假設，認為變量X與Y相關，且高度正相關。

(7) 協同性檢驗（Kendall相關檢驗）

選取1999年到2008年運動員獲世界冠軍男女人數資料，試用Kendall檢驗方法檢驗二者之間是否存在相關關系？

計算結果如下：

Kendall's rank correlation tau

data:  x and y

T = 39, p-value = 0.002213

alternative hypothesis: true tau is not equal to 0

sample estimates:

      tau

0.7333333

結果顯示P值=0.002681 < 0.05，是以拒絕原假設，認為變量X與Y相關，且高度正相關。

(8) 多總體比較的秩和檢驗

本案例資料來源于百度文庫。

某汽車駕駛員記錄了使用5種不同牌子的汽油每5加侖行駛的距離，資料如下：

品牌1	37.5	31.3	33.8	32.5
品牌2	36.3	32.5	36.3	35.0
品牌3	40	40	43.8	46.3
品牌4	36.3	42.5	40	41.3
品牌5	40	32.5	38.8	33.8

這些資料是否說明這5種牌子的汽油每加侖平均行駛距離全相等?

計算結果如下：

Kruskal-Wallis rank sum test

data:  x

Kruskal-Wallis chi-squared = 11.996, df = 4, p-value = 0.01738

結果顯示P值=0.017 < 0.05，故拒絕原假設，認為這5種牌子的汽油每加侖平均行駛距離不全相等。

(9) 多總體比較的卡方檢驗

本案例資料來自課後習題

在研究教師的疲勞對教學效果的影響時，有人做了如下實驗，讓一名教師在一天中對一個班中的三個組教授同一門課程，這三個組是随機劃分的，最後的考試成績為：

學生(成績)	1	2	3
8點上課	72	84	75
10點上課	77	67	79
14點上課	58	78	80

試問每組間的學生的學習有無顯著差異，取顯著性水準為0.05。

計算結果如下：

Pearson's Chi-squared test

data:  x

X-squared = 4.4667, df = 4, p-value

= 0.3465

結果顯示P值=0.3465 > 0.05，故接受原假設，認為每組間的學生的學習無顯著差異。

實驗項目三

 (1)制作資料經驗分布函數,使用分布拟合方法解決總體類型的檢驗問題

本例給出15名學生的體重資料（機關：kg）

75.0	64.0	47.4	66.9	62.2
62.2	58.7	63.5	66.6	64.0
57.0	69.0	56.9	50.0	72.0

繪制出15名學生體重的經驗分布圖和相應的正态分布圖：

繪制出的經驗分布圖和正态分布曲線：

将學生體重進行排序得：

47.4 50.0 56.9 57.0 58.7 62.2 62.2 63.5 64.0 64.0 66.6 66.9 69.0 72.0 75.0

求得經驗分布函數為：

(2) 機率密度圖像，使用分布拟合方法解決總體類型的檢驗問題

繪制出直方圖和密度估計曲線和正态分布的機率密度曲線：

通過上圖，可以明顯看出密度估計曲線和正态分布的機率密度曲線還是有一定的差别的。

實驗項目四

通過最小二乘與權函數結合使用的方法解決非線性回歸問題

本案例資料來自課後習題

1.一隻紅鈴蟲的産卵數與溫度有關。下表是産卵數Y與溫度X的一組資料，試研究Y與X的回歸關系。

編号	1	2	3	4	5	6	7
溫度x	21	23	25	27	29	32	25
産卵數y	7	11	21	24	66	115	325

繪制出散點圖：

X	Y	線性拟合	殘差	權函數估計	混合估計
21	7	50.91986	-43.919861	-43.919861	7
23	11	63.06620	-52.066202	-52.066202	11
25	21	75.21254	-54.212544	97.787456	173
27	24	87.35889	-63.358885	-63.358885	24
29	66	99.50523	-33.505226	-33.505226	66
32	115	117.72474	-2.724739	-2.724739	115
25	325	75.21254	249.787456	97.787456	173

解得回歸直線方程為：

附錄

#非參數統計

#1

data1 <- read.csv("C:\\Users\\Administrator\\Desktop\\運動員.csv",header=T)

cham_num <- data1[,3]

cham_item <- data1[,2]

list("人數排序"=sort(cham_num),"五數總括"=fivenum(cham_num),"分位數"=quantile(cham_num))

x <- cham_item

list("均值"=mean(x),"截尾均值"=mean(cham_num,trim = 0.05))

 qj<-function(p,alpha){

   n=length(x);

   for(i in 1:n){

     s1=ppois(i-1,n*p);

     s2=ppois(i,n*p); #s1+dpois(i,n*p)

     if(s1<=(1-alpha)/2&&s2>(1-alpha)/2) break

     }

   for(j in n:1){

     s3=1-ppois(j,n*p);

     s4=1-ppois(j-1,n*p); #s3+dpois(j,n*p)

     if(s3<=(1-alpha)/2&&s4>(1-alpha)/2) break

     }

   dp<-paste("(",i,",",j,")")

   ci<-paste("[",sort(x)[i],",",sort(x)[j],"]")

   list("資料位置"=dp,"區間估計"=ci)

   }

 qj(0.5,0.9)

x<-cham_item

y<-cham_num

list("樣本均值之差"=mean(y)-mean(x),"樣本中位數之差"=median(y)-median(x))

 HL<-function(x,y,alpha){

   n1<-length(x);

   n2<-length(y);

   d<-n1*n2/2-sqrt(n1*n2*(n1+n2+1)/12)*qnorm(1-alpha/2);

   nn<-c();

   for(j in 1:n2){

     nn<-c(nn,y[j]-x)

     }

   d1=floor(d+1);d2=n1*n2-d1;

   wz<-paste(c(d1,d2+1))

   ci<-paste("[",sort(nn)[d1],",",sort(nn)[d2+1],"]")

   list("Yi-Xi資料"=nn,sort(nn),中位數=median(nn),上下界位置=wz,"95%置信區間"=ci)

   }

HL(x,y,alpha=0.05)

#2

#中位數檢驗

sign.test<-function(x, m0, alpha=0.05, alter="two.sided"){

   p<-list( )

   n<-length(x)

   sign<-as.numeric(x>=m0)

   s<-sum(sign)

   result<-binom.test(s, n, p=0.5, alternative=alter,

                        conf.level=alpha)

   p$p.value=result$p.value

   p

   }

 sign.test(y,155)

#兩總體比較的Wilcoxon秩和檢驗

data2 <- read.csv("C:\\Users\\Administrator\\Desktop\\運動員男女獲得世界冠軍.csv",header=T)

cham_man <- data2[,2]

cham_woman <- data2[,3]

x<-cham_man

y<-cham_woman

z<-append(x,y)

rxy<-rank(z);rxy

ry<-rank(z)[-c(1:length(x))];ry

jrs<-sum(ry);jrs

r1<-function(m,n,alpha){round(m*(m+n+1)/2-qnorm(1-alpha/2)*sqrt(m*n*(m+n+1)/12))}

r2<-function(m,n,alpha){round(m*(m+n+1)/2+qnorm(1-alpha/2)*sqrt(m*n*(m+n+1)/12))}

list("r1"=r1(10,13,0.05),"r2"=r2(10,13,0.05),"是否落入拒絕域"=l<-r1(10,13,0.05)<jrs&jrs<r2(10,13,0.05)) #是否落入接受域

#兩總體比較的K-S檢驗

ks.test(x,y)

#卡方獨立性檢驗

yesbelt <- c(12813,647,359,42)

nobelt <- c(65963,4000,2642,303)

chisq.test(rbind(yesbelt,nobelt))

#Fisher列聯表檢驗

x <- c(12813,65963,1048,6945);dim(x) <- c(2,2)

fisher.test(x)

#Spearman秩相關檢驗方法

x<-data2[1:10,2]

y<-data2[1:10,3]

cor.test(x,y,method = "spearman")

#協同性檢驗（Kendall相關檢驗）

x<-data2[1:10,2]

y<-data2[1:10,3]

cor.test(x,y,method = "kendall")

#多總體比較的秩和檢驗

x <- list(

  a=c(37.5,31.3,33.8,32.5),

  b=c(36.3,32.5,36.3,35.0),

  c=c(40,40,43.8,46.3),

  d=c(36.3,42.5,40,41.3),

  e=c(40,32.5,38.8,33.8)

)

kruskal.test(x)

#多總體比較的卡方檢驗

x <- c(72,77,58,84,67,78,75,79,80)

dim(x) <- c(3,3)

x

chisq.test(x,correct = FALSE)

#制作資料經驗分布函數,使用分布拟合方法解決總體類型的檢驗問題

w <- c(75,64,47.4,66.9,62.2,62.2,58.7,63.5,66.6,64,57,69,56.9,50,72)

plot(ecdf(w),verticals = TRUE,do.p=FALSE)

x <- 44:78

lines(x,pnorm(x,mean(w),sd(w)))

sort(w)

# 機率密度圖像，使用分布拟合方法解決總體類型的檢驗問題

w

hist(w,freq = F,main="機率密度直方圖")

lines(density(w),col="blue")

x<-44:78

lines(x,dnorm(x,mean(w),sd(w)),col="red")

#試研究Y與X的回歸關系

x<-c(21,23,25,27,29,32,25)

y<-c(7,11,21,24,66,115,325)

plot(x,y)

ML<-function(x,y,bw){

   a<-lm(y~x) #線性拟合但隻傳回拟合系數用unclass()具體檢視

   b<-a$fitted.values #提取拟合結果清單

   c<-a$residuals #拟合殘差結果清單

   d<-matrix(c(x,c),nc=2)

   e<-c()

   for(i in 1:length(x)){

     f<-which((d[i,1]-bw)<d[,1]&d[,1]<(d[i,1]+bw))

     e<-c(e,mean(d[f,2]))

    }

   data.frame(X=x,Y=y,線性拟合=b,殘差=c,權函數估計=e,混合估計=b+e)

}

ML(x,y,0.5)

a<-{mean(y)*mean(x^2)-mean(x*y)*mean(x)}/{mean(x^2)-mean(x)^2};a

b<-{mean(x*y)-mean(x)*mean(y)}/{mean(x^2)-mean(x)^2};b