非參數(shù)統(tǒng)計(jì)習(xí)題R程序

1、P37.例 2.1 build.price-c36,32,31,25,28,36,40,32,41,26,35,35,32,87,33,35 ;build.price 60,88,88,87,60,73,60,97,91,60,83,87,81,90;length scores# 輸入向量求長(zhǎng)度ss-cscores-80;ss t-0 t1-0 fori in 1:lengthss histbuild.price,freq=FALSE#直方圖linesdensitybuild.price,col=red#連線# 方法一： m-meanbuild.price;m#均值if ssi0 t-t+1#

2、求小于 80 的個(gè)數(shù)D-varbuild.price#方差else t1-t1+1求大于 80 的個(gè)數(shù)SD t;t1 1 13 1 15 binom.testsumscores80,lengthscores,0.75 t=m-37/SD/sqrtlengthbuild.price;t#t統(tǒng)計(jì)量運(yùn)算檢驗(yàn)統(tǒng)計(jì)量 t= 1 -0.1412332 # 方法二： t.testbuild.price-37#課本第 38 頁(yè)p-value = 0.0014360.01 例 2.2 Cox-Staut 趨勢(shì)存在性檢驗(yàn)P47 binom.testsumbuild.price37,lengthbuild.pric

3、e,0.5# 課本 40 頁(yè) 例 2.3 P-2*1-pnorm1.96,0,1;P 例 2.6 year-1971:2022;year lengthyear 1 0.04999579 P1-2*1-pnorm0.7906,0,1;P1 rain 例 2.4 p t1=0 p1-2*pnorm-0.9487,0,1;p1 for i in 1:5 1 0.3427732 if rainiraini+5 t1-t1+1 t1 例 2.5（ P45）k-0:t1-1 scores-c95,89,68,90,88,60,81,67,60,60,60,63,60,92, sumdbinomk,5,0.

4、5# =0.1875 y-6/25;y# =0.1875 #2 該地區(qū)前 32 年降雨量是否變化？例 2.9 t=0 shuju39-data.frameread.table for i in 1:16 SHUJU39.txt,header=TRUE;shuju39 if rainiraini+16 t-t+1 attachshuju39 sum.a=0 t sum.b=0 k1-0:mint,16-t-1 sum.c=0 sumdbinomk1,16,0.5# =0.0002593994 for i in 1:lengthid pbinommaxk1,16,0.5#= 0.000259399

5、4 y1-1+16/216;y1#=0.0002593994 plotyear,rain ablinev=1971+2022/2,col=2 if pinzhongi=A sum.a-sum.a+chanliangi else if pinzhongi=B sum.b-sum.b+chanliangi else fuhao-sum.c-sum.c+chanliangi linesyear,rain anovalmrainyear sum.a;sum.b;sum.c ma-sum.a/4 mb-sum.b/4 mc-sum.c/4 ma;mb;mc 隨機(jī)游程檢驗(yàn)（P50）例 2.8 client

6、-cF,M,M,M,M,M,F,M, M,F,M,M,M,M,F,M,F, fuhao-repa,12;fuhao M,M,M,F,F,F,M,M,M;client for i in 1:lengthid n0 n1-sumclient=M;n1 n0-n-n1;n0 t1-0 for i in 1:lengthclient-1 if clienti=clienti+1 t1-t1 else t1-t1+1 fuhaoi0 fuhaoi0 fuhaoi-+ else fuhaoi- fuhao R-t1+1;R#=12 #find rejection region（不寫(xiě)）# 利用上題編程解決檢

7、驗(yàn)的隨機(jī)性rl-1+2*n1*n0/n1+n0*1-1.96/sqrtn1+n0;rl n-lengthfuhao;n ru-2*n1*n0/n1+n0*1+1.96/sqrtn1+n0;ru#=15.33476 （課本為 ru=17 ）n1-sumfuhao=+;n1 n0-n-n1;n0 t1-0 for i in 1:lengthfuhao-1 if fuhaoi=fuhaoi+1 t1-t1 p else t1-t1+1 N=19 #sample size of expected distribution of W+ R-t1+1;R y-dwilxonfunN;y #find rej

8、ection region rl-1+2*n1*n0/n1+n0*1-1.96/sqrtn1+n0;rl ru-2*n1*n0/n1+n0*1+1.96/sqrtn1+n0;ru # 運(yùn)算 PW+=x 中的 x 取值的 R 參考程序！dwilxonfun=functionN a=c1,1 #when n=1 frequency of W+=1 or o 例 2.10 （P52）libraryquadprog# 不存在叫n=1 pp=NULL #distribute of all size from 2 to N aa=NULL #frequency of all size from 2 to

9、N for i in 2:N t=crep0,i,a a=ca,rep0,i+t p=a/2i #density of Wilcox distribut when size=N a quadprog 這個(gè)名字的程輯包libraryzoo# 不存在叫zoo這個(gè)名字的程輯包librarytseries# 不存在叫tseries 這個(gè)名字的程輯包 run1=factorc1,1,1,0,rep1,7,0,1,1,0,0,rep1,6,0,r ep1,4, 0,rep1,5,rep0,4,rep1,13；run1 y=factorrun1 runs.testy# 錯(cuò)誤 : 沒(méi)有 runs.test 這

10、個(gè)函數(shù) N=19 #sample size of expected distribution of W+ Wilcoxon 符號(hào)秩檢驗(yàn)y-dwilxonfunN;lengthy-1 W+ 在零假設(shè)下的精確分布W+ 分布密度函數(shù),histy,freq=FALSE # 下面的函數(shù)dwilxonfun用來(lái)運(yùn)算linesdensityy,col=red 即 PW+=x 的一個(gè)參考程序！例 2.12 （P59）dwilxonfun=functionN a=c1,1 #when n=1 frequency of W+=1 or o ceo-c310,350,370,377,389,400,415,425,

11、440,295, 325,296,250,340,298,365,375,360,385;lengthceo n=1 # 方法一pp=NULL #distribute of all size from 2 to N wilcox.testceo-320 aa=NULL #frequency of all size from 2 to N # 方法二 for i in 2:N ceo.num320;ceo.num t=crep0,i,a n=lengthceo a=ca,rep0,i+t binom.testceo.num,n,0.5 p=a/2i #density of Wilcox dist

12、ribut when size=N 例 2.13P61 a-c62,70,74,75,77,80,83,85,88 walsh-NULL for i in 1:lengtha-1 for j in i+1:lengtha walsh-cwalsh,ai+aj/2 walsh=cwalsh,a NW=lengthwalsh;NW printcol col1-crstuconffl,rstuconffl+1,conffl+2;col1 min 例 2.14 “stu-c82,53,70,73,103,71,69, 80,54,38,87,91,62,75,65,77;stu alpha=0.05

13、n=lengthstu;n medianwalsh conf=pbinomn,n,0.5-pbinom0,n,0.5;conf fork in 1:n 2.5 單組數(shù)據(jù)的位置參數(shù)置信區(qū)間估量P61 conf=pbinomn-k,n,0.5-pbinomk,n,0.5 例 2.14 if conf1-alphaloc=k-1;break stu-c82,53,70,73,103,71,69, 80,54,38,87,91,62,75,65,77;stu printloc alpha=0.05 （剩余的例題參考程序在課本）rstu-sortstu;rstu 3.6 正態(tài)記分檢驗(yàn)conff-NULL

14、;conff 例 2.18 n=lengthstu;n fori in 1:n-1 baby11-alphaconff-cconff,i,j,conf 例 2.18 qiuzhi-functionx n=lengthx conff a=rep2,n lengthconff for i in 1:n min-103-38;min ai=sumx=xi c-seq1,lengthconff-1,3;c fori in c a col-crstuconffi,rstuconffi+1,conffi+2 min1-rstuconffi+1-rstuconffi fuhao-functionx,y if

15、 min1minmin-min1;l-i n=lengthx sgn=rep2,n pvaluey sgni=1 例 2.20 else if xi=y shuju2-data.framechongshu=c0:6, sgni=0 zhushu=c10,24,10,4,1,0,1;shuju2 else attachshuju2 sgni=-1 n=sumzhushu;n lamda-sumchongshu*zhushu/n;lamda sgn p-dpoischongshu,lamda;p n*p n1-lengthbaby x.squ=sumzhushu2/n*p-n;x.squ baby

16、zhi=qiuzhibaby # 方法一q=n1+1+babyzhi/2*n1+2 value-qchisq1-0.05,lengthzhushu-1;value babysgn-fuhaobaby,34 # 方法二babysgn=signbaby1-34;babysgn s=qnormq,0,1 pvalue-1-pchisqx.squ,lengthzhushu-1;pvalue W-ts%*%babysgn;W 例 2.21 sd-sums*babysgn2;sd T=W/sd;T shuju3-c36,36,37,38,40,42,43,43,44,45,48,48, 50,50,51,

17、52,53,54,54,56,57,57,57,58,58,58,58, 58,59,60,61,61,61,62,62,63,63,65,66,68,68,70, 2.7 分布的一樣性檢驗(yàn) 例 2.19 shuju1-data.framemonth=c1:6, 73,73,75;shuju3 n=lengthshuju3 n0=sumshuju330 & shuju340 & shuju3=50;n2 n-sumcustomers;n expect50 & shuju360 & shuju370 & shuju380;n6 運(yùn)算 45 位同學(xué)value-qchisq1-0.05,length

18、customers-1;value nn-cn0,n1,n2,n3,n4,n5,n6;nn # 方法二體重分類的頻數(shù)！shuju3.mean=meanshuju3;shuju3.mean shuju3.var=varshuju3;shuju3.var #md.xy-quantilexy,0.25 #利用 p 分位數(shù)的檢shuju3.sd=sdshuju3;shuju3.sd 驗(yàn)e0=pnorm30,shuju3.mean,shuju3.sd tmd.xy e1=pnorm40,shuju3.mean,shuju3.sd- lx-lengthx pnorm30,shuju3.mean,shuju

19、3.sd ly-lengthy e2=pnorm50,shuju3.mean,shuju3.sd- lxy-lx+ly pnorm40,shuju3.mean,shuju3.sd e3=pnorm60,shuju3.mean,shuju3.sd- Amd.xy pnorm50,shuju3.mean,shuju3.sd if alt=greater e4=pnorm70,shuju3.mean,shuju3.sd- w-1-phyperA,lx,ly,t pnorm60,shuju3.mean,shuju3.sd e5=pnorm80,shuju3.mean,shuju3.sd- else i

20、f alt=less pnorm70,shuju3.mean,shuju3.sd w-phyperA,lx,ly,t e6=1-pnorm80,shuju3.mean,shuju3.sd e=ce0,e1,e2,e3,e4,e5,e6;e ee=n*ce0,e1,e2,e3,e4,e5,e6;ee x.squ=sumnn2/ee-n;x.squ conting.table=matrixcA,lx-A,lx,t-A,ly-t-A,ly,t,lxy-t,lxy,3,3 -cX,Y,X+Y MXY,MXY,TOTAL # 方法一 dimnamesconting.tab

21、, value-qchisq1-0.05,lengthee-1;value listcontingency.table=conting.table,p.vlue=w # 方法二例 3.2 pvalue-1-pchisqx.squ,lengthee-1;pvalue X-c698,688,675,656,655,648,640,639,620 Y-c780,754,740,712,693,680,621 例 2.22 healthy-c87,77,92,68,80,78,84,77,81,80,80,77,92,86, 76,80,81,75,77

22、,72,81,90,84,86,80,68,77,87,76,77,7 8,92, 75,80,78;healthy ks.testhealthy,pnorm,80,6 第三章# 方法一：BM.testX,Y,less # 方法二：XY-cX,Y md.xy-medianXY tmd.xy lx-lengthX ly-lengthY lxy-lx+ly Amd.xy # 沒(méi)有修正時(shí)的情形 pvalue1-pnormA,lx*t/lx+ly, #Brown_Mood中位數(shù)#Brown-Mood中位數(shù)檢驗(yàn)程序BM.test-functionx,y,alt xy-cx,y md.xy-medianx

23、y #利用中位數(shù)的檢驗(yàn)sqrtlx*ly*t*lx+ly-t/lx+ly3;pvalue1 # 修正時(shí)的情形 pvalue2-pnormA,lx*t/lx+ly-0.5, Mx-My 的 R 參考程序：x1-c140,147,153,160,165,170,171,193 x2-c130,135,138,144,148,155,168 sqrtlx*ly*t*lx+ly-t/lx+ly3;pvalue2 n1-lengthx1 n2-lengthx2 th.hat-medianx2-medianx1 3.2 、Wilcoxon-Mann-Whitney秩和檢驗(yàn)# 求兩樣本分別的秩和的程序. Q

24、iuzhi-functionx,y n1-lengthy B=10000 Tboot=crep0,1000 #vector of length Bootstrap yyyy,1 wm 例 3.3 weight.low=c134,146,104,119,124,161, replacement from x1 xx2=samplex2,5,T #sample of size n2 with replacement from x2 Tbooti=medianxx2-medianxx1 th-medianTboot;th se=sdTboot 107,83,113,129,97,123 m=leng

25、thweight.low weight.high=c70,118,101,85,112,132,94 n=lengthweight.high Normal.conf=cth+qnorm0.025,0,1*se,th-qnorm0.025,0,1*se;Normal.conf Percentile.conf=c2*th-quantileTboot,0.975,2*th-quantile # 方法一：Tboot,0.025;Percentile.conf wy-Qiuzhiweight.low,weight.high#wy=50 Provotal.conf=cquantileTboot,0.025

26、,quantileTbo ot,0.975;Provotal.conf wxy-wy-n*n+1/2;wxy#=22 mean-m*n/2 var-m*n*m+n+1/12 pvalue-1-2*pnormwxy,mean-0.5,var;pvalue th.hat 3.3 、Mood 方差檢驗(yàn) qiuzhi-functionx,y xy-cx,y # 方法二zhi-NULL for i in 1:lengthx wilcox.testweight.high,weight.low zhi=xy 例 3.4 zhi 3.4 、Moses 方差檢驗(yàn) 引例：qiuzhi-functionx,y x1

27、-c48,56,59,61,84,87,91,95 xy-cx,y x2-c2,22,49,78,85,89,93,97 zhi-NULL zhi_x1=qiuzhix1,x2;zhi_x1 for i in 1:lengthx #zhi_x2=qiuzhix2,x1;zhi_x2 zhi=xy #var_x1=varx1;var_x1 #var_x2=varx2;var_x2 zhi m=lengthx1;m n=lengthx2;n 例 3.6 mean_R=m+n+1/2;mean_R x1-c8.2,10.7,7.5,14.6,6.3,9.2,11.9, mean1=m*m+n+1*m

28、+n-1/12;mean1 5.6,12.8,5.2,4.9,13.5 var1=m*n*m+n+1*m+n+2*m+n-2/180;var1 m1=lengthx1;m1 M1=sumzhi_x1-mean_R2;M1 x2-c4.7,6.3,5.2,6.8,5.6,4.2, p_value=2*pnormM1,mean1-0.5,sqrtvar1 6.0,7.4,8.1,6.5 隨機(jī)分組p_value m2=lengthx2;m2 例 3.5 A-matrixx1,ncol=3;A#a1=samplex1,3,F X-c4.5,6.5,7,10,12 xx2=NULL Y-c6,7.2,8

29、,9,9.8 fori in 1:m1 zhi_X=qiuzhiX,Y;zhi_X ifsuma1=x1i=0 xx2=cxx2,x1i m=lengthX;m n=lengthY;n a2=samplexx2,3,F mean_R=m+n+1/2;mean_R mean2=m*m+n+1*m+n-1/12;mean2 xx3=NULL var2=m*n*m+n+1*m+n+2*m+n-2/180;var2 fori in 1:m1-3 M2=sumzhi_X-mean_R2;M2 ifsuma2=xx2i=0 xx3=cxx3,x1i # 方法一：查附表9 # 方法二：a3=samplexx

30、3,3,F p_value=2*1-pnormM2,mean2-0.5,sqrtvar2 x11=sumA1,-meanx12;x11 p_value # 方法三x12=sumA2,-meanx12;x12 x13=sumA3,-meanx12;x13 Z=1/sqrtvar2*M2-mean2+0.5;Z x14=sumA4,-meanx12;x14 SSA-cx11,x12,x13,x14;SSA x-factorlever;x B-matrixx21:9,ncol=3;B xy-data.framey,x y11=sumB1,-meanx22;y11 attachxy y12=sumB2

31、,-meanx22;y12 y13=sumB3,-meanx22;y13 aovformula=yx,data=xy aov.xy-aovformula=yx,data=xy SSB-cy11,y12,y13;SSB summaryaov.xy zhi_SSA=qiuzhiSSA,SSB;zhi_SSA zhi_SSB=qiuzhiSSB,SSA;zhi_SSB # 方法二：運(yùn)算總的平方和S=sumzhi_SSA;S x1-c1.4,1.9,2.0,1.5 TM=S-4*4+1/2;TM x2-c2.0,2.4,1.8,2.2 # 方法一 查附表 4 x3-c2.6,2.8,2.5,2.1 拒

32、絕域 C=TMW0.025,m1,m2或者yW0.975,m1,m2 y.mean-meany;y.mean 其中 W0.975,m1,m2=m1*m2-W0.025,m1,m2. ssT-sumy-y.mean2;ssT # 方法二（ Wilcoxon 秩和檢驗(yàn)）x1.mean-meanx1 wilcox.testSSA,SSB x2.mean-meanx2 # 方法二（ Mann-Whitney秩和檢驗(yàn)）x3.mean-meanx3 m=lengthSSA;m n=lengthSSB;n sse-sumsumx1-x1.mean2, sumx2-x2.mean2,sumx3-x3.mean

33、2;sse # 運(yùn)算誤差平方和mean_AB=m*n/2;mean_AB var_AB=m*n*m+n+1/12;var_AB sst-ssT-sse;sst #運(yùn)算組間平方和運(yùn)算方差分析的Fp_value=1-F-sst/2/sse/lengthy-3;F #pnormS,mean_AB,sqrtvar_AB;p_value 檢驗(yàn)統(tǒng)計(jì)量第四章# 臨界值的運(yùn)算4.1 、試驗(yàn)設(shè)計(jì)和方差分析的基本概念回憶value-qf0.95,2,lengthy-3;value #R 軟件中單因素方差分析的函數(shù) # 運(yùn)算 p-value 值例 4.1 p.value-1-pf8,2,lengthy-3;p.va

34、lue # 方法一：表 4.5 *Analysis of Variance Model * y-c2.0,1.4,2.0,2.8,2.4,1.9,1.8,2.5,2.0,1.5,2.1,2.2;y lever-cB,A,C,C,B,A,B,C,A,A,C,B xueye-c8.4,9.4,9.8,12.2, 10.8,15.2,9.8,14.4,8.6,9.8,10.2,9.8, 8.8,9.8,8.9,12.0,8.4,9.2,8.5,9.5;xueye sst1-sumxueye-meanxueye2;sst1 a=matrixxueye,ncol=5;a quzu-applya,2,su

35、m;quzu chuli-applya,1,sum;chuli k=5 H=12/n*n+1*sumazhi2/lengtha+sumbzhi 2/lengthb +sumczhi2/lengthc+sumdzhi2/lengthd-3*n+1 b=4 方法一： value=qchisq1-0.05,3;value ssb=1/4*sumquzu2-sumquzu2/k*b;ssb sst=1/5*sumchuli2-sumchuli2/k*b;sst sse=sst1-ssb-sst;sse 方法二： pvalue=1-pchisqH,3;pvalue mean=cmeana,meanb,m

36、eanc,meand mssb=ssb/k-1;mssb # 兩兩比較的程序msst=sst/b-1;msst msse=sse/k*b-k-b+1;msse F1=mssb/msse;F1 F2=msst/msse;F2 value1=qf1-0.05,k-1,k*b-k-b+1 value2=qf1-0.05,b-1,k*b-k-b+1 bjiao=functionazhi,bzhi,czhi,dzhi n=lengthcazhi,bzhi,czhi,dzhi av=sumazhi/lengthazhi bv=sumbzhi/lengthbzhi se=sqrtn*n+1/12*1/len

37、gthazhi+1/lengthbzhi d=absav-bv dab=d/se huizong=cd,se,dab,qnorm1-0.05,0,1 例 4.3 qiuzhi-functionw,x,y,z xy-cw,x,y,z zhi-NULL for i in 1:lengthw zhi=xy zhi a-c80,203,236,252,284,368,457,393 b-c133,180,100,160 c-c156,295,320,448,465,481,279 huizong bjiaoazhi,bzhi,czhi,dzhi bjiaoczhi,dzhi,azhi,bzhi 4.3

38、 、Jonckheere-Terpstra檢驗(yàn)例 4.5 x=c125,136,116,101,105,109 y=c122,114,131,120,119,127 z=c128,142,128,134,135,131,140,129 xm=meanx;xm ym=meany;ym zm=meanz;zm d-c194,214,272,330,386,475 g=crep1,6,rep2,6,rep3,8 azhi=qiuzhia,b,c,d;azhi tapplycx,y,z,g,median bzhi=qiuzhib,a,c,d;bzhi JT.testdata=tcx,y,z,class=g czhi=qiuzhic,a,b,d;czhi Wij-functionx,y dzhi=qiuzhid,a,b,c;dzhi n1=length