北京科技大學(xué)概率論與數(shù)理統(tǒng)計上機報告3

1、概率論與數(shù)理統(tǒng)計第三次上機報告專信息與計算科學(xué)業(yè)班級學(xué)生姓名指導(dǎo)教師張志岡U完成時間信計1502 （35組）呂瑞杰 陳炎睿 何芝芝2020年1月16日UniversiiyScience cind TechnDbogy BeijingMatlab概率論與數(shù)理統(tǒng)計上機練習(xí)（3）五、假設(shè)檢驗【例】（離散型分布檢驗）某工廠近五年發(fā)生了 63起事故，按星期幾可以分為 9 10 11 8 13 12，問該廠發(fā)生的事故數(shù)是有與星期幾有關(guān)？clear allmi=9 10 11 8 13 12; %周一到周六的事故數(shù) n=sum(mi); %總的事故數(shù)r=0; %總體中沒有未知參數(shù) k=length(mi);

2、 % 天數(shù) pii=1/6; %事故的概率 kai2=0;kai2=sum(mi -n*pii).A2)./(n*pii);% k2 統(tǒng)計量的值alpha1=0.05; %顯著性水平 alpha2=0.01; %顯著性水平 alpha3=0.001; %顯著性水平 la1=chi2i nv(1 -alpha1,k-r-1); % kai2 la2=chi2i nv(1 -alpha2,k-r-1); % kai2 la3=chi2i nv(1 -alpha3,k-r-1); % kai2 pz=1-chi2cdf(kai2,k -r-1);%右側(cè)概率 if kai2la2xzx=*;分布的累計

3、概率，即臨界值 分布的累計概率，即臨界值 分布的累計概率，即臨界值elseif kai2la1xzx=*;elsexzx=-;endx=0:0.1:la3; y=chi2pdf(x,k -r-1); plot(x,y);x1=kai2:0.1:la3; y1=chi2pdf(x1,k -r-1); hold onif kai2=60);se6=le ngth(fi nd(se=60);%及格人數(shù)sy7=le ngth(fi nd(sy=80);se7=le ngth(fi nd(se=80);%優(yōu)秀人數(shù)sy8=sy6/sy1;se8=se6/se1;% 及格率sy9=sy7/sy1;se9=s

4、e7/se1;% 優(yōu)秀率fprintf（t人數(shù)t 平均分t 最小值t最大值 t極差tt標(biāo)準(zhǔn)差tt及格人數(shù)及格率t優(yōu)秀人數(shù)優(yōu)秀率n）;fprintf（n);fprintf(數(shù)學(xué) 4dt%10.4ft%4dtt%4dtt%4dtt%10.4ft%4dt%10.4ft%4dt%10.4fn,sy1,sy4,sy2min,sy2max,sy3,sy5,sy6,sy8,sy7,sy9)fprintf(信計 %4dt%10.4ft%4dtt%4dtt%4dtt%10.4ft%4dt%10.4ft%4dt%10.4fn,se1,se4,se2min,se2max,se3,se5,sy6,sy8,se7,se

5、9)fprin tf(n);%方法一fprintf(檢驗數(shù)學(xué)和信計的方差是否相等n);h1,p1,varci1,stats1=vartest2(sy,se,alpha,both);if(h仁=0)disp(結(jié)果：方差相等);elsedisp(結(jié)果：方差不相等);endfprin tf(n);% %方法二% F=sy5A2/se5A2;%統(tǒng)計量F,滿足F分布% alpha=0.05;%取顯著水平為 0.05% Fla1=finv(alpha/2,sy1 -1,se1-1);Fla2=finv(1 -alpha/2,sy1-1,se1-1);% 求 F 的臨界值% if (FFla1 & FFla

6、2)% MM=數(shù)學(xué)分析1和數(shù)學(xué)分析2的方差無顯著差異；% else% MM=數(shù)學(xué)分析1和數(shù)學(xué)分析2的方差有顯著差異；% end% fprintf(檢驗數(shù)學(xué)分析1和數(shù)學(xué)分析2的方差是否相等n);% fprintf(統(tǒng)計量F的值ttt顯著性水平tt臨界值ttttt檢驗結(jié)果n);% fprintf( %.4ftttt%.4fttt%.4fttt%15sn,F,alpha,Fla1,MM);% fprin tf(nn);%方法一fprintf(檢驗數(shù)學(xué)和信計的平均分是否相等n);h2,p2,muci2,stats2=ttest2(sy,se,alpha,both);if(h2=0)disp(結(jié)果：平均

7、分相等);elsedisp（結(jié)果：平均分不相等）;endfprin tf(n);%方法二% %方法三% sw=(sy1 -1)*sy5A2+(se1 -1)*se5A2)/(sy1+se1 -2);% T=(sy4-se4)/sw/sqrt(1/sy1+1/se1);% 統(tǒng)計量 T,滿 T 分布% Tla1=tinv(alpha/2,sy1+se1 -2);Tla2=tinv(1 -alpha/2,sy1+se1-2);% 求出 T 的臨界值% if (abs(T)Ua)disp(優(yōu)秀率無顯著差異);elsedisp(優(yōu)秀率有顯著差異);endP=(sy6+se6)/(sy1+se1);U=(

8、sy8 -se8)/sqrt(sy8+se8)*p*(1 -p);Ua=norminv(1 -alpha/2);if(abs(U)Ua)disp(及格率無顯著差異);elsedisp(及格率有顯著差異);endh3,p3,kstat3,critval3=lillietest(sy,alpha); if(h3=1)disp(數(shù)學(xué)不是正態(tài)分布)elsedisp(數(shù)學(xué)是正態(tài)分布) endh4,p4,kstat4,critval4=lillietest(se,alpha); if(h4=1)disp(信計不是正態(tài)分布)elsedisp(信計是正態(tài)分布)end%hist(sy)% 直方圖%h5,p5,s

9、tats5=chi2gof(sy)% 可以檢驗分布%cdf=sy ,no rmcdf(sy,sy4,sy5)%h5,p5,ksstat,cv5=kstest(sy,cdf)% a=0:1:100;% a=a;% CDF=a,cdf(a,sy4,sy5);% h = kstest(sy,CDF,0.05);S=6565688174766869827774667372776062816668766074809069606368676962606060676077676060607172606661866460608973744340619569706266636278746050766265847

10、06983734371707371697460616070747848936461797153606052636061656278606560858589696660;h,p,jbstat,critval=jbtest(S,alpha);if(h=0)disp(服從正態(tài)分布);else disp(不服從正態(tài)分布);endsavg=mea n( S);svar=var(S);x=20:130;y=no rmpdf(x,savg,sqrt(svar);d=5;a=20:d:130;pdf=hist(S,a)./le ngth(S)./d;plot(x,y,r);hold onscatter(a,p

11、df,filled);hold off輸出： Ix3_1_lrj_41521335人數(shù)平均分 最小值 最大值極差標(biāo)準(zhǔn)差及格人數(shù)及格率優(yōu)秀人數(shù)優(yōu)秀率數(shù)學(xué) 5470.666740955510.0885530.981590.16679.2313530.9815信計 6065.516743894650.0833檢驗數(shù)學(xué)和信計的方差是否相等結(jié)果：方差相等檢驗數(shù)學(xué)和信計的平均分是否相等結(jié)果：平均分不相等優(yōu)秀率有顯著差異及格率有顯著差異數(shù)學(xué)不是正態(tài)分布Warning: P is less than the smallest tabulated value, returning 0.001. In lilli

12、etest (line 206)In Ix3_1_lrj_41521335 (line 99)信計不是正態(tài)分布服從正態(tài)分布逼 Figure 1一 XFile Edit View jnstrl Tools Desktop Wi ndcw M*lp a jT 謠，s a六、方差分析【例1】(單因素方差分析)考慮溫度對某化工產(chǎn)品得率的影響，選擇五種不同溫度進行試驗，每一溫度各做三次試驗。 方法一：自編程序clear allX=90,92,88;97,93,92;96,96,93;84,83,88;84,86,82;a=5;ni=3,3,3,3,3;%每個因素的樣本數(shù)n=sum( ni);%樣本總數(shù)%

13、T=sum(sum(X); %先求每列的和，再求總和%求所有樣本的和 T,平方和，及ST,SA,SET=0; ST=0; SA=0; SE=0;for i=1:aTi=0;forj=1: ni (i)T=T+X(i,j);Ti=Ti+X(i,j);ST=ST+X(i,j)A2; endSA=SA+T22/ni(i);%總偏差平方和%效應(yīng)平方和%誤差平方和endST=ST-TA2/n;SA=SA-TA2/n;SE=ST-SA;F=(SA/(a-1)/(SE/( n-a); alpha1=0.05;%顯著性水平la1=fi nv(1-alpha1,a-1, n-a); alpha2=0.01;%顯

14、著性水平la2=fi nv(1-alpha2,a-1, n-a);%由卩分布的累積概率，求臨界值,PFla2%由卩分布的累積概率，求臨界值，%計算F比值做為臨界點的右側(cè)概率p=1-PXFPFla1xzx=elsexzx=endfprintf( fprintf( fprintf( fprintf(來源tt 平方和tt 自由度tt 均方和ttF效應(yīng) Att%.2ftt%4dtt%.2ftt%.2ftt%.4fn誤差 tt%.2ftt%4dtt%.2ftttt%4sn總和 tt%.2ftt%4dtt臨界值比tt 顯著性n);,SA,a-1,SA/(a-1),F,p);,SE ,n-a,SE/( n-

15、a),xzx);=%.2f(%.2f),%.2f(%.2f)n,ST ,n-1,la1,alpha1,la2,alpha2);fprintf( nn);運行結(jié)果為：來源平方和自由度均方和F比顯著性效應(yīng)A303.60475.9015.18*誤差50.00105.000.000299總和353.6014臨界值=3.48(0.05),5.99(0.01)方法二：調(diào)用 matlab工具anoval(X)，其中矩陣 X表示X的轉(zhuǎn)置，即該函數(shù)每一列為一個因素。 運行結(jié)果為【例7.2(沒有交互作用的多因素方差分析)一火箭使用了四種燃料，三種推進器，作射程試驗。X=58.2,56.2,65.3;49.1,54

16、.1,51.6;60.1,70.9,39.2;75.8,58.2,48.7;方法一：自編程序，運行結(jié)果為58.200056.200065.300049.100054.100051.600060.100070.900039.200075.800058.200048.7000來源平方和自由度均方和F比顯著性效應(yīng)A157.59352.530.4306X (0.738747)效應(yīng)B223.852111.920.9174X (0.449118)誤差731.986122.00總和1113.4211臨界值=4.76(0.05),5.14(0.05)方法二：調(diào)用 matlab工具anova2(X)Matlab

17、程序?qū)崿F(xiàn):X=60,60,63,63,40,69,65,60,72,67,62,78,82,90,69,60,72,76,78,93,69,68,95,71,83,60; 73,73,60,74,77,71,85,70,89,60,61,77,62,68,60,70,66,84,74,69,61,60,86,73,69,74;50,81,67,65,77,71,76,62,89,65,65,62,62,60,78,81,66,70,80,53,69,66,61,48,156,69;6i8,60,74,60,62,43,61,60,60,64,70,74,65,7379,60,43,76,66,

18、63,60,60,68,60,60,60;ano val(X)a=4;ni=26,26,26,26;%每個因素的樣本數(shù)n=sum( ni);%樣本總數(shù)%T=sum(sum(X);%先求每列的和，再求總和%求所有樣本的和 T,平方和，及 ST,SA,SET=0; ST=0; SA=0; SE=0;for i=1:a nTi=0;for j=1: ni (i)T=T+X(i,j);Ti=Ti+X(i,j);ST=ST+X(i,j)A2;endSA=SA+TiA2/ni(i);endST=ST-TA2/n; %總偏差平方和SA=SA -TA2/n;%效應(yīng)平方和SE=ST-SA; %誤差平方和F=(S

19、A/(a -1)/(SE/(n-a);% F 比alpha1=0.05; %顯著性水平Ia1=finv(1 -alpha1,a-1,n-a); % 由 F 分布的累積概率，求臨界值，PFla=1 -alphaalpha2=0.01; %顯著性水平la2=finv(1 -alpha2,a-1,n-a); % 由 F 分布的累積概率，求臨界值，PFla=1 -alphap=1-fcdf(F,a-1,n-a);%計算F比值做為臨界點的右側(cè)概率p=1 -PXla2xzx=*;elseif Fla1xzx=*;elsexzx=-;endfprintf(ttt對四個班的數(shù)學(xué)分析一nn)fprintf(來源

20、tt平方和ttt自由度ttt均方和tttF比tt顯著性n);fprintf(效應(yīng) Att%.2ftttt%4dtt%.2ftt%.2ftt%.4fn,SA,a-1,SA/(a-1),F,p);fprintf(誤差 tt%.2ftt%4dtt%.2fttttt%4sn,SE,n -a,SE/(n-a),xzx);fprintf(總和 tt%.2ftt%4dtt 臨界值=%.2f(%.2f),%.2f(%.2f)n,ST,n -1,la1,alpha1,la2,alpha2); fprin tf(nn);x=6060747275 6481607673767873956373818485958563

21、486674605074404032462554396996988477919365617270606073607673897580677287898794979467858490859493628078646780787810096939297978289888090909290879485939785697679698186856061827460716572718170666778766083806351847885837884928393999784908596696075776772726861637274858595998491879589719493819487968377777

22、585908660637768687985739288758288797385848788909660787880768889749991879092927796929091949671728787858088859998969594977084878381838989918290919196606977896073826172818575928877808475868591627177736680876888828289818960606669665073706071807292906664708671789384100848790959474878185918481698181828878

23、926164606370817360656573616674861009887969088736477818683926960727375787874748685959091717375838280787471787776769350626770687874819094869190916790918588859065524356656276779993909095967145647263836676739388927490628577778787898994939391929565523985688260656140857674686237455364557162856391779188604

24、939685071497892809193789281999888959199668783938285917071828378879080998992969496536675746962786981789191759466778582849579616044616061664834636441666066503672607163696363726660716149327253786060606379607081607167898368828586808060826552364553464852687189696186836045458339756074608083889081608775899

25、390896250307364767143606568627060616869727679776064678670777960867875758383646067714673607069677964778074697666658068657980878380857392898486849679748382748590606062837387734360678380907076779689969196666070726174726377898877928760486856687171606163666560636876827581847760516367397369606782837886866

26、0647379728186678789878788977473698690959564626482799185；X=x：a=7;ni=114,114,114,114,114,114,114;% 每個因素的樣本數(shù)n=sum( ni); %樣本總數(shù)%T=sum(sum(X);%先求每列的和，再求總和%求所有樣本的和 T,平方和，及 ST,SA,SET=0; ST=0; SA=0; SE=0;for i=1:aTi=0;for j=1: ni (i)T=T+X(i,j);Ti=Ti+X(i,j);ST=ST+X(i,j)A2;endSA=SA+TiA2/ni(i);endST=ST-TA2/n; %

27、總偏差平方和SA=SA -TA2/n;%效應(yīng)平方和SE=ST-SA; %誤差平方和F=(SA/(a -1)/(SE/(n-a);% F 比alpha1=0.05; %顯著性水平Ia1=finv(1 -alpha1,a-1,n-a); % 由 F 分布的累積概率，求臨界值，PFla=1 -alphaalpha2=0.01; %顯著性水平la2=finv(1 -alpha2,a-1,n-a); % 由 F 分布的累積概率，求臨界值，PFla=1 -alphap=1-fcdf(F,a-1,n-a);%計算F比值做為臨界點的右側(cè)概率p=1 -PXla2xzx=*;elseif Fla1xzx=*;el

28、sexzx=-;endfprintf(ttt 全體數(shù)學(xué)課程 nn)fprintf(來源tt平方和ttt自由度tt均方和tttF比tt顯著性n);fprintf(效應(yīng) Att%.2ftt%4dtt%.2ftt%.2ftt%.4fn,SA,a-1,SA/(a-1),F,p);fprintf(誤差 tt%.2ftt%4dtt%.2fttttt%4sn,SE,n-a,SE/(n-a),xzx);fprintf(總和 tt%.2ftt%4dtt 臨界值=%.2f(%.2f),%.2f(%.2f)n,ST,n -1,la1,alpha1,la2,alpha2); fprin tf(nn);輸出：方法一：對

29、四個班的數(shù)學(xué)分析一來源平方和自由度均方和F比顯著性效應(yīng)A906.263302.093.130.0289誤差9636.2710096.36*總和10542.53103臨界值=2.70(0.05),3.98(0.01)全體數(shù)學(xué)課程來源平方和自由度均方和F比顯著性效應(yīng)A15804.9362634.1515.340.0000誤差135824.95791171.71*總和151629.87797臨界值=2.11(0.05),2.82(0.01)方法* Figure 1: One-way ANOVA XFile drt view nsert Tools 衛(wèi)專號 ktop Wi ndow Help四ANOV

30、A TableScure tssdf監(jiān)FFrcbFCaluanj;屈.3-S73. LB. 0289Erzor9E36. 3iaa9fi 3 83Total：32石 Figure 2-口XFile Edit View nsert Tools DesktopWi ndew Help1已il給 Q七、回歸分析【例1】（一元線性回歸）以三口之家為單位，某食品在某年中平均月消費量（kg）與其價格（元/kg）之間的關(guān)系。方法一：自編程序clear allx=2,2,2.4,2.7,3,3,3.5,3.6,3.8,4,4.5,5;y=3,3.6,2.8,2.8,2.3,2.9,1.9,2.1,1.9,1.

31、3,1.5,1;運行結(jié)果為回歸直線方程為y=4.8256+ -0.7799 x來源平方和自由度均方和F比顯著性回歸R6.040016.040090.2608*誤差0.6692100.06690.000003總和6.709211臨界值=4.9646(0.05),10.0443(0.01)【練習(xí)3.3回歸分析(1) 對全體學(xué)生的“基礎(chǔ)外語一”和“基礎(chǔ)外語二”進行回歸分析；(2) 對全體學(xué)生的“數(shù)學(xué)分析一”和“體育一”進行回歸分析。“基礎(chǔ)外語一”和“基礎(chǔ)外語二”回歸直線方程為y=(11.8263)+(0.8605)x來源平方和自由度均方和F比顯著性回歸 R 6613.480016613.480024

32、9.1623誤差2972.800711226.5429p=0.000000總和9586.2807113 臨界值=3.9258(0.05),6.8667(0.01)Figure 1匚1 迂|反File Edit View Ins ert呂 曲疲 t 町 Tiaiiw Jfelp S A畋匡1H H回歸直線方程為 y=(92.6417)+( -0.0551)x來源平方和自由度均方和F比顯著性回歸 R34.5457134.54570.5277總和7366.7368113臨界值=3.9258(0.05),6.8667(0.01)Matlab程序?qū)崿F(xiàn)：基礎(chǔ)外語I和基礎(chǔ)外語IIx=7783 78 87 7

33、5 81 83 67 74 83 81 82 74 74 94 65 81 77 80 87 8376 85 87 76 79 83 75 89 87 69 95 80 65 69 84 74 70 87 607089837160556671736369706752737172826964626865697167537270646867866473647877787873787968687268768983668175738772728160698460706076806574867885637258;y=8076 728977878466798980817980937587888984817

34、78692808581758994719070576788797887616988856863566575776872737060796169837366706672727670607381726173868075768481708266757869677078758982656364789475748262628563686271776969828687697756;n=len gth(x);X=ones(n,1) x;alpha=0.05;b,b in t,r,ri nt,stats=regress(y,X,alpha);ahat=b(1);bhat=b(2);hold ontitle($

35、haty=hat11.8263+hat 0.8605x$,i nterpreter,latex,fo ntsize,20) yhat=b(1)+b(2).*x; scatter(x,y,filled); plot(x,yhat,k);for i=1: n;a=x(i) x(i); c=yhat(i) y(i); plot(a,c,r-);end sxx=sum(x.A2)-(sum(x)A2/n;t=bhat*sqrt(sxx)/sqrt(stats(4);p=(1-tcdf(abs(t),n-2)*2;se1= stats(4);sr=sum(yhatmea n(y)42);st=sum(y

36、-mea n( y)42);se=sum(yhat-y).A2);disp(stats(4);for i=1: n;y1(i)=yhat(i)-(-tinv(0.025,n-2)*sqrt(1+1/n+(x(i) -mean(x).A2/sxx)*sqrt(se/(n-2);y2(i)=yhat(i)+( -tinv(0.025,n-2)*sqrt(1+1/n+(x(i) -mean(x).A2/sxx)*sqrt(se/(n-2); end yleft=y1;yright=y2;Y=yleft yright;plot(x,yleft,g,x,yright,g)hold offf=sr/(se

37、/( n-2);lamda1=fi nv(0.95,1, n-2);lamda2=fi nv(0.99,1, n-2);p=1-fcdf(f,2-1,n-2);%計算F比值做為臨界點的右側(cè)概率p=1-PXlamda1)a=;endif(flamda2&flamda1)a= endif(flamda2) a=-;end fprintf(來源tt平方和tt自由度tt均方和ttF比tt顯著性n)fprintf(回歸 tt%.4ft%.ftt%.4ft%.4ftt%sn,sr,1,sr,f,a)fprintf(誤差tt%.4ft%.ftt%.4ftp=%.8fn,se,n -2,se/(n-2),p)

38、fprintf(總和 tt%.4ft%.f臨界值=%.4f(0.05),%.4f(0.01)n,st,n-1,lamda1,lamda2)輸出:來源平方和自由度均方和F比顯著性回歸6613.480016613.4800 249.1623*誤差2972.800711226.5429p=0.00000000總和9586.2807113臨界值=3.9258(0.05),6.8667(0.01)凱 Figure 1一 XFile Edit View nstrt Tools Desktop Wi ndcw N*lp、勺=11.8263 + 0.8605110 -100 -50 -4011111115055606570750085 盹 95數(shù)學(xué)分析I和體育Ix=656568817476686982777466737277606281666876607480906960636867696260606067607767606