《數(shù)據(jù)科學(xué)統(tǒng)計計算》課后習(xí)題答案chap4

統(tǒng)計計算第四章作業(yè)

L設(shè)X1,X2,...,X“和匕,心,…,4分別是來自總體N3，b'2)和N(〃2,b2)(b未知)的兩個獨

立樣本，檢驗原假設(shè)和備擇假設(shè)分別為：/：=〃2,"1：*〃2.取n=10(),7?7=50,Px=

0,〃2=2,。=1,用數(shù)值計算方法估計上述檢驗問題的功效。

解答：統(tǒng)計量/y為：

Txv——1-?tin+iTi—2),

+ZH-1)

其中

sw=(Qx+。>.)/("+m-2),

I2

QX=2(x，一幻之，。產(chǎn)x(匕-n.

Z=11=1

數(shù)值計算方法估計功效如下

I

mul=0

mu2=2

n=100

m=50

sigma=1

K=10000

res<-c()

fbr(iin1:K){

X<-rnorm(n,mu1,sigma)

Y<-rnorm(m,mu2,signia)

Sw<-(sum((X-mean(X))A2)+sum((Y-mean(Y))A2))/(n+m-2)

T<-(mean(X)-mean(Y))/sqrt(Sw*(l/n+1/m))

res[i]v-as.numeric(abs(T)>=qt(0.975,n+m-2))#取絕對值進行判斷

}

result<-mean(res)#數(shù)值模擬估計的統(tǒng)計功效

criti<-qt(0.975,n+m-2);

power<-2-pt((criti一(mul-mu2)/sqrt(Sw*(l/n+l/m))),m+n-2)-pt((criti+(mul-mu2)/sqrt(Sw*(l/n+1Zm))),m+n

-2)

#理論統(tǒng)計功效

c(result,power)

此例子中，理論統(tǒng)計功效和數(shù)值模擬的的結(jié)果均為1。

2.35位健康男性在未進食前的血糖濃度如下，分別用卡方檢驗和單樣本K-S檢驗方法來檢

驗這組數(shù)據(jù)是否來自均值〃=80、標(biāo)準差。=6的正態(tài)分布總體，用數(shù)值計算方法進行分析。

數(shù)據(jù)：6868727575767677777777777778787880808080808081818184848686

878790929292

解答：(1)卡方檢驗：需要將樣本數(shù)據(jù)進行類別劃分。對于本例35個樣本數(shù)據(jù)，可以劃分

為4個類別，對于正態(tài)分布來說就是劃分為4個區(qū)間，如果數(shù)據(jù)落入這4個區(qū)間的概率值是相

等的，也就是每個區(qū)間的概率為1/4=0.25o查標(biāo)準正態(tài)分布表，可得到標(biāo)準正態(tài)分布下劃分4個

區(qū)間的3個臨界值點。相應(yīng)的可以得出均值為〃=80、標(biāo)準差夕=6的三個臨界值所對應(yīng)的X

值(區(qū)間界限)，如下表：

概率值對應(yīng)z值對應(yīng)X值

0.75-0.6775.98

0.500.0080.00

0.250.6784.02

0.00正無窮正無窮

根據(jù)計算的區(qū)間界限，將樣本數(shù)據(jù)的每一個數(shù)據(jù)都分配到對應(yīng)的區(qū)間中，統(tǒng)計實際頻數(shù)。最后由

實際頻數(shù)與期望頻數(shù)的差值計算卡方統(tǒng)計量，計算過程如下表：

劃分區(qū)間Oi實際頻數(shù)Ei期望頻數(shù)(Oi-Ei)2/Ei

75.98以下58.751.607143

75.98-80.00118.750.5785714

80.00-84.02118.750.5785714

84.02以上88.750.06428571

合計35352.828572

根據(jù)卡方檢驗內(nèi)容可知該統(tǒng)計量服從自由度為的卡方分布。本例中，劃分了4個區(qū)間k=4,自

由度為4-1=3。查卡方分布表(單側(cè)顯著水平0.05,自由度為3)得卡方值為。因為計算卡方統(tǒng)

計量為2.829<7.815,落在接受域，所以接受原假設(shè)，拒絕備擇假設(shè)，即觀測數(shù)據(jù)可以認為服從

正態(tài)分布。

(2)單樣本K-S檢驗的算法如下

n=35

data2=c(68,68,72,75,75,76,76,77,77,77,77,77,77,78,78,78,80,80,80,80,80,80,81,81,

81,84,84,86,86,87,87,90,92,92,92)

data2<-sort(data2)

Dsplus<-max(abs(c(l:n)/n-pnonn(data2,mean=80,sd=6)))

Dminus<-max(abs((c(1:n)-1)/n-pnorm(data2,mean=80>sd=6)))

2

statI<—max(D_splus,D_minus)

c(statl,0.22424)

由于0.148小于臨界值0.224,不拒絕原假設(shè)，則認為樣本服從均值為80,標(biāo)準差為6的正態(tài)分

布，與卡方檢驗的結(jié)果一致。

3.設(shè)Xi,X2,Xn和Ki,Y2,分別是來自總體-—+1)和U(02,■+1)的兩個獨立樣

本，檢驗原假設(shè)和備擇假設(shè)分別為：“°：仇=。2,4：01大出。取n=30,m=50A=0,02=0.5分

別用曼-惠特尼U檢驗和兩樣本游程檢驗方法分析上述假設(shè)檢驗問題，并采用數(shù)值計算方法估計

功效。

解答：(1)曼-惠特尼U檢驗數(shù)值分析如下

set.seed(422)

n=30

m=50

theta1=0

theta2=0.5

K=1000

resu=NULL

fbr(iin1:K){

X_dala<-runif(n,theta1,theta1+1)

Y_data<-runif(m,thcta2,thcta2+1)

Data=c(X_data,Y_data)

Rank=rank(Data)

W_x=sum(Rank[l:n])

W_y=sum(Rank[(n+1):(n+m)])

Wxy=(W_x-n*(n+l)/2)/(m*n)

stat=(W_xy-1/2)/sqrt((m+n-1)/(12*m*n))

resu[i]=as.numeric(abs(stat)>1,96)

)

mean(resu)

數(shù)值模擬結(jié)果是1,曼-惠特尼U檢驗具有良好的檢驗功效。

(2)兩樣本游程檢驗數(shù)值分析如下

###同分布檢驗-Wald-Wolfowitz游程檢驗

set.seed(422)

n=30

m=50

N=m+n

theta1=0

theta2=0.5

K=1000

3

resl=NULL

fbr(iin1:K){

X_data<-runif(n,theta1,thetal+1)

Y_data<-runif(m,theta2,theta2+1)

#原假設(shè)thetal=theta2

Data=c(X_data,Y_data)

Runs=as.numeric(order(Data)<=n)

W=1

for(jin1:(length(Runs)-1)){

if(Runsfj]!=Runs[j+l])

W=W+1

)

E=2*m*n/N;

V=2*m*n*(2*m*n-N)/(NA2*(N-1))

#resu[i]=(W-E)/sqrt(V)

resl[i]=as.numeric(abs((W-E)/sqrt(V))>=1.96)

)

mcan(resl)

i

兩樣本游程檢驗的數(shù)值檢驗功效為0.995o

4.設(shè)Xi,X2,X”和匕,力,…,Yni分別是來自總體X和Y的兩個獨立樣本,X和Y的分布函

數(shù)分別為F(x)和G(y)，檢驗的原假設(shè)和備擇假設(shè)分別為：

為:E(x)=G(y),“i:F(x)，G(y)

取〃?=50,”=30,X?N((),1),/?N((),2),兩樣本K-S檢驗構(gòu)造檢驗統(tǒng)計量，并采用數(shù)值計

算方法估計檢驗統(tǒng)計量的功效。

解答：兩樣本K-S檢驗數(shù)值算法如下

###同分布檢驗-K-S檢驗(柯氏檢驗)

set.seed(422)

m=50

n=30

N=m+n

statl=NULL

resl=NULL

res2=NULL

fbr(iin1:1000){

#創(chuàng)建數(shù)據(jù)

X_data<-rnonn(n,0,1)

Y_data<-rnorm(m,0,sqrt(2))

#計算每個Xi,Yi的經(jīng)驗分布

F_nx<-((matrix(1,n,l)%*%X_data)<=(X_data%*%matrix(l,1,n)));F_nx[F_nx]<-1;Fnx<-rowMeans(

F_nx)

4

F_ny<-((matrix(1,m,1)%*%X_data)<=(Y_data%*%matrix(1,1,n)));F_ny[F_ny]<-l;F_ny<-rowMeans(

F_ny)

Grnx<-((matrix(1,n,1)%*%Y_data)<=(X_data%*%matrix(1,1,m)));G_mx[G_mx]<-1;G_mx<-

rowMcans(G_mx)

G_my<-((matrix(1,m,1)%*%Y_data)<=(Y_data%*%matrix(1,1,m)));G_my[G_my]<-1;G_my<-

rowMeans(G_my)

D_x=max(abs(F_nx-G_mx))

D_y=max(abs(F_ny-G_my))

statl=max(D_x,D_y)

resl[i]=as.numeric(statl>=1.36*sqrt(N/(m*n)))#確定臨界值

index=scq(1,10000,1)

pval=2*sum((-1)A(index-1)*exp(-2*m*n*indexA2*stat1A2/N))

res2[i]=as.numcric(p_val<0.05)

#判斷P值是否小于犯一類錯誤的概率

)

c(mean(res1),mean(res2))

5.分類變量X和Y分別有r和c種可能性，分別記為X==1,2,…,r)以及V=/(/=

l,2,...,c).將同時滿足X=i(i=l,2,...,r)和丫=j(J=l,2,...,c)的樣本觀測頻數(shù)記為固定

r,c,利用列聯(lián)分析方法檢驗X和Y是否獨立，并用數(shù)值分析方法計算檢驗的功效。

解答:

Y1Y2???Yc合計

XInilnl2nlcnl.

X2n21n22n2cn2.

:

Xrnrlnr2nrcnr.

合計

對書上P66的具體問題進行復(fù)現(xiàn)與代碼注釋：

取樣本量〃=100,假設(shè)數(shù)據(jù)(為,匕),(x2,%),…,(X”,乙)服從二元正態(tài)分布，x和y的均值

均為0,方差都為1,相關(guān)系數(shù)為0.6+0.04*p,P為1到10的整數(shù)。

將兩組數(shù)據(jù)劃分為k=8個區(qū)間，計算檢驗的功效。

###獨立檢驗列聯(lián)表檢驗

library(MASS)

set.seed(422)

efficacy4<-c()

fbr(pin1:10){

k=8

n=100

Sigma=matrix(c(1,0.6+0.04*p,0.6+0.04*p,l),2,2)

resl=NULL

fbr(iin1:1000){

5

data=mvmorm(n,rep(0,2),Sigma)#dala自動生成兩列數(shù)據(jù)

X_data=data[,1]#將第一列賦值給X

Y_data=data[,2]#將第二列賦值給Y

Sx=sort(Xdata)#排序后方便后續(xù)縱向比較

Sy=sort(Ydata)

#根據(jù)每組數(shù)據(jù)極值劃分區(qū)間

interx=seq(min(Xdata),max(Xdata),by=(max(X_data)-min(Xdata))/k)

inter_y=seq(min(Y_data),max(Y_data),by=(max(Y_data)-min(Y_data))/k)

#構(gòu)造X的區(qū)間

left_inter_x=rcp(l,n)%*%t(inter_x[l:k])

right_inter_x=rep(1,n)%*%t(inter_x[2:(k+1)])

#構(gòu)