統(tǒng)計建模與R軟件 第四講 (2015)

1、關(guān)于統(tǒng)計量的誘導(dǎo)關(guān)系：關(guān)于統(tǒng)計量的誘導(dǎo)關(guān)系：22211() ,( ,)1niiiSXXXNn (0,1)N(0,1)iN22( ),( )nm2(0,1),( )Nn22S2St221( )niin22(1)2N/( ,)/nF n mm2F( )/t nn2Nt2N) 1(1)(ntSnn) 1(222nSnn兩個正態(tài)母體誘導(dǎo)的統(tǒng)計量：兩個正態(tài)母體誘導(dǎo)的統(tǒng)計量：222(,)iN 211(,)iN 兩個完全不同的正態(tài)分布母體誘導(dǎo)F分布22(,)iN 21(,)iN 211 12212222 211(1)(1,1)(1)nnn SnF nnn Sn具有相同方差的正態(tài)分布母體誘導(dǎo)t分布122211

2、222122nnn Sn SSnn121212()()(2)11t nnSnn主要內(nèi)容主要內(nèi)容4.1 矩法4.2 極大似然估計4.3 估計量的優(yōu)良性準(zhǔn)則4.4 區(qū)間估計思想：用樣本矩去估計總體矩，總體矩與總體的參數(shù)有關(guān)，從而得到總體參數(shù)的估計。設(shè)總體X的分布函數(shù)F(x;1 m )中有m個未知參數(shù), 假設(shè)總體的m階原點矩存在,n個樣本x1 xn ，令總體的k階原點矩等于樣本的k階原點矩，即4.1 4.1 矩法矩法解此方程組得到則稱 為參數(shù)k的矩法估計量。m1,.,k一階，二階矩法估計參數(shù)：一階，二階矩法估計參數(shù)：更一般的提法為:利用樣本的數(shù)字特征作為總體的數(shù)字特征的估計.例如:無論總體服從什么分

3、布，其均值和方差分別為:解得均值與方差的矩法點估計:2211()1niiSXXn20(1)01kMkMppp 設(shè)總體服從二項分布B(k; p);k, p為未知參數(shù)。X1,x2,xn是總體X的一個樣本，求參數(shù)k,p的矩估計 。,k pM1是總體均值(一階原點矩)M2是總體方差(二階中心矩)1Mkp2(1)kppM 解得：kR實現(xiàn)實現(xiàn):(1)# N=20,p=0.7, 試驗次數(shù)n=100 x m11 13.84 m21 4.8544# 由解析計算給定結(jié)果：N=m12/(m1-m2); N #1 21.31695 p=(m1-m2)/m1; p # 1 0.6492486kR實現(xiàn)實現(xiàn):(2)1221

4、(1)fkpMffkppMmoment_fun-function(p) f-c(p1*p2-M1, p1*p2-p1*p22-M2) J-matrix(c(p2, p1, p2-p22, p1-2*p1*p2),nrow=2, byrow=T) list(f=f, J=J)1122(1)2ffpkkpJffppkkpkp牛頓法：Newtons-function (fun, x, ep=1e-5, it_max=100) index-0; k-1 while (k=it_max) x1 - x; obj - fun(x); x - x - solve(obj$J, obj$f); norm -

5、sqrt(x-x1) %*% (x-x1) if (normep) index-1; break k-k+1 obj - fun(x); list(root=x, it=k, index=index, FunVal= obj$f)#fun是列表,返回函數(shù)表達(dá)式和函數(shù)的Jacobi矩陣;x是迭代初值#初始化#把初值記下來#牛頓法:x1=x0-J-1f#index是示性指標(biāo)，如果等于1表示牛頓法解存在，否則沒有解#函數(shù)返回一個列表:根,迭代次數(shù),示性指標(biāo),函數(shù)值主函數(shù)：x-rbinom(100, 20, 0.7); n-length(x)M1-mean(x); M2-(n-1)/n*var(x)s

6、ource(moment_fun.R); source(Newtons.R)p-c(10,0.5); Newtons(moment_fun, p)f,JNewtons-function (fun, x, ep=1e-5, it_max=100) index-0; k-1 while (k=it_max) x1 - x; obj - fun(x); x - x - solve(obj$J, obj$f); norm - sqrt(x-x1) %*% (x-x1) if (normep) index-1; break k-k+1 obj - fun(x); list(root=x, it=k, i

7、ndex=index, FunVal= obj$f)K0,p0$root1 20.9158983 0.6564385$it1 520(1)01kMkMppp極大似然法極大似然法定義1：設(shè)總體X的概率密度函數(shù)或分布律為是未知參數(shù)， 為來自總體X的樣本，稱( , ),f x 12,nXXX121( ; )( ;,)( , )nniiLxLx xxf x ( ();)sup ( ;)LXXLX為的似然函數(shù)(likelihood function).定義2：設(shè)總體X的概率密度函數(shù)或分布律為是未知參數(shù)， 為來自總體X的樣本，為的似然函數(shù),若： 是一個統(tǒng)計量，且滿足：( , ),f x 12,nXXX(

8、; )Lx為12X)= (X ,X ,X )n（則稱 為的極大似然估計.1.似然函數(shù)關(guān)于似然函數(shù)關(guān)于連續(xù)連續(xù)極值條件，得：( ;)0iLX 似然方程。獨立同分布的樣本，似然函數(shù)具有連乘的形式ln ( ;)0iLX 例子：正態(tài)分布例子：正態(tài)分布對數(shù)似然方程：2211()niiXXn11niiXxn# multiroot()函數(shù)計算# e1=mu, e2=sigma, x=樣本model mean(x)1 0.1273094 sum(x-mean(x)2)/101 1.267102$root1 0.2480794 0.907706422221ln ( ,; )ln(2)()22inLxx 121(

9、)0iLFx222241()022iLnFx 2.似然函數(shù)關(guān)于似然函數(shù)關(guān)于有間斷點有間斷點設(shè)總體X服從區(qū)間a; b的均勻分布，x=x1; ; xn為來自總體的一組樣本，用極大似然估計法估計參數(shù)a; b。似然函數(shù)為L(a; b,x)不是a; b的連續(xù)函數(shù),其似然方程為:不能求解從極大似然估計的定義出發(fā)來求L(a; b,x)的最大值,要使L達(dá)到最大，那么b-a應(yīng)該盡可能的小，但是a不能大于min(x),b不能小于max(x),因此a; b的極大似然估計為:3.是離散參數(shù)空間是離散參數(shù)空間一般地：在魚塘釣出r條魚，做上記號，然后再釣出s條，發(fā)現(xiàn)有x條有標(biāo)記第二次釣出的魚的條數(shù)x服從超幾何分布：(;

10、)()(1; )L N xg NL Nx()s xxN srsNCCP XxC似然函數(shù)為L(N;x)=P(X=x)似然函數(shù)的比為：()()()Ns NrN Nsrx 22()()Nrs NrsNrs NxNrsNx將數(shù)字帶入上式得池塘中魚的總數(shù)為：500*1000/72=6944例子：例子：在魚塘撈出500條魚，做上記號，然后再撈出1000條，發(fā)現(xiàn)有72條有標(biāo)記，試估計魚塘所有的魚有多少？4.在解似然方差時無法得到解析解，采用數(shù)值方法在解似然方差時無法得到解析解，采用數(shù)值方法設(shè)總體X服從Cauchy分布，其概率密度函數(shù)為：21111 ()nniix21( ; ),1 () f xxx 其中為未

11、知參數(shù).X1,X2,Xn是總體X的樣本,求極大似然估計.Cauchy分布的似然函數(shù)為：1( ; )( , )niiLxf x2ln ( ; )lnln(1 () )iLxnx 求導(dǎo)2101 ()niiixx求對數(shù)似然方程的解析解是困難的，考慮使用數(shù)值方法。1.使用uniroot函數(shù):# 參數(shù)為1的cauchy分布x=rcauchy(100,1)f-function(p) sum(x-p)/(1+(x-p)2)out out$root1 0.9020655$f.root1 1.800204e-072.使用OPTMIZE()函數(shù)x=rcauchy(100,1)loglike optimize(lo

12、glike,c(0,5)2# ln ( ; )lnln(1 () )iLxnxminimum =0.9021objective =254.4463exitflag =1#的近似解#-lnL(,x)的近似值$minimum1 1.03418$objective1 239.4626matlab解#-lnL=min,則lnL=max, #optimize只能求最小，最大問題轉(zhuǎn)化為負(fù)的最小問題關(guān)于二項分布的極大似然估計：關(guān)于二項分布的極大似然估計：( , ;)(1)iiiiXXn XnXL n p XCpp ( , ;)()(1)maxiiiiXmnXXnXL n p XCpp ()(1)iiiXXn

13、 Xinp xXCppmatlab輸出的極大似然估計數(shù)值解：x = 20.0000 0.7065fval = 210.2846%matlabfunction f=fg(sita)x=load(abc.txt);s=0;for i=1:100 s1=log(nchoosek(fix(sita(1),x(i); s2=log(sita(2)*x(i); s3=log(1-sita(2)*(sita(1)-x(i); s=s+s1+s2+s3;endf=-s;%matlab主函數(shù)：x0=21,0.5A=0,1;0,-1;-1,0b=1;0;-20 x,fval = fmincon(fg,x0,A,b

14、)矩法估計值：$root1 20.9158983 0.6564385$it1 5ln ( , ;)(ln)iXnL n p XC ()ln(1)imnXplniXpR實現(xiàn)：實現(xiàn)：obj=function(n)x-rbinom(100, 20, 0.7); m-length(x)f=-sum(log(choose(n1%/%1),x)-(log(n2)*sum(x)+log(1-n2)*(m*n1-sum(x)sita0=c(20,0.5)#初值constrOptim(sita0, obj, NULL, ui=rbind(c(0,-1),c(0,1),c(1,0), ci=c(-1,0,-20)

15、ln ( , ;)(ln)iXnL n p XC ()ln(1)imnXplniXpR輸出結(jié)果：$par1 22.0340214 0.6179089$value1 209.5277matlab輸出的極大似然估計數(shù)值解：x = 20.0000 0.7065fval = 210.2846結(jié)果對比區(qū)間估計：區(qū)間估計：設(shè)總體X的分布函數(shù)F(x,)含有未知參數(shù)，對于給定值(0 1),若由樣本X1,X2,Xn確定的兩個統(tǒng)計量， 和 滿足:112(,)nXXX212(,)nXXX112212( (,)(,)1nnPXXXXXX 則稱隨機區(qū)間 是參數(shù)的置信度為1- 的置信區(qū)間。12 ( ,) 置信下限置信上限

16、置信度越短越好3.1一個正態(tài)總體的情況一個正態(tài)總體的情況3.1.1均值的區(qū)間估計:2 已知時：(0,1)/XNn/2/2,XXnZZn/21/XPZn 參數(shù)的置信度為1-的雙側(cè)置信區(qū)間2 未知時：(1)/Xt nSn/2/2(1)(1),SStntnXXnn/21/XPtSn 參數(shù)的置信度為1-的雙側(cè)置信區(qū)間interval_estimate1-function(x,sigma=-1,alpha=0.05) n-length(x); xb=0) tmp-sigma/sqrt(n)*qnorm(1-alpha/2); df-n else tmp-sd(x)/sqrt(n)*qt(1-alpha/

17、2,n-1); df-n-1 data.frame(mean=xb, df=df, a=xb-tmp, b=xb+tmp)#默認(rèn)未知#函數(shù)返回一個數(shù)據(jù)框例子：例子： 4.14 某工廠生產(chǎn)的零件長度X被認(rèn)為服從N(,0.04),現(xiàn)從該產(chǎn)品中隨機抽取6個，其長度的測量值如下(單位：mm):試求該零件長度的置信系數(shù)為0.95的區(qū)間估計。15.115.214.814.915.114.6source(interval_estimate1.R)x=c(14.6,15.1,14.9,14.8,15.2,15.1)interval_estimate1(x,sigma=0.2) mean df a b 14.9

18、5 6 14.78997 15.11003 置信區(qū)間t.test(x): One Sample t-testdata: x t = 162.1555, df = 5, p-value = 1.692e-10alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: 14.71300 15.18700 sample estimates:mean of x 14.95假設(shè)：H1拒絕H1 (接受H0)的概率幾乎所有的統(tǒng)計軟件P-value都是這個意思當(dāng)已知時：3.1.2方差方差 的區(qū)間估計的區(qū)

19、間估計22221()( )niiXn22n22/22( )1nPn 參數(shù) 的置信度為1-的雙側(cè)置信區(qū)間22222/21/2,( )( )nnnn當(dāng)未知時：222(1)(1)nSn222/21/22(1)(1)1),(1)nSnSnn22/22(1)(1)1nSPn 參數(shù) 的置信度為1-的雙側(cè)置信區(qū)間2interval_var1-function(x,mu=Inf,alpha=0.05) n-length(x) if (muInf) S2 - sum(x-mu)2)/n; df - n else S2 - var(x); df - n-1 a-df*S2/qchisq(1-alpha/2,df)

20、 b-df*S2/qchisq(alpha/2,df) data.frame(var=S2, df=df, a=a, b=b)#默認(rèn)未知，未知標(biāo)志=inf#已知時， muInf#未知時， mu=Inf例例4.16：用區(qū)間估計方法估計例4.15的測量誤差2，分別對均值已知(=10)和均值未知兩種情況進(jìn)行討論。#輸入數(shù)據(jù)，調(diào)用編好的程序x=c(10.1,10,9.8,10.5,9.7,10.1,9.9,10.2,10.3,9.9)interval_var1(x,mu=10) var df a b0.055 10 0.02685130 0.1693885interval_var1(x)var df

21、a b0.05833333 9 0.02759851 0.19441644.3.2兩個正態(tài)總體的情況兩個正態(tài)總體的情況interval_estimate2-function(x, y, sigma=c(-1,-1),var.equal=FALSE, alpha=0.05) n1-length(x); n2-length(y) xb-mean(x); yb=0) #均值差1- 2的區(qū)間估計(置信度為1-) tmp-qnorm(1-alpha/2)*sqrt(sigma12/n1+sigma22/n2) df-n1+n2 else if (var.equal = TRUE) Sw-(n1-1)*v

22、ar(x)+(n2-1)*var(y)/(n1+n2-2) tmp-sqrt(Sw*(1/n1+1/n2)*qt(1-alpha/2,n1+n2-2) df-n1+n2-2 else S1-var(x); S2-var(y) nu-(S1/n1+S2/n2)2/(S12/n12/(n1-1)+S22/n22/(n2-1) tmp-qt(1-alpha/2, nu)*sqrt(S1/n1+S2/n2) df t.test(x,y) #兩個樣本方差不同 Welch Two Sample t-testdata: x and y t = 0.7551, df = 24.467, p-value = 0

23、.4574alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -1.594229 3.436546 sample estimates:mean of x mean of y 500.8576 499.9365 t.test(x,y,var.equal=TRUE)2.配對數(shù)據(jù)情形下均值差的區(qū)間估計配對數(shù)據(jù)情形下均值差的區(qū)間估計編號12345678910X11.315.015.013.512.810.011.012.013.012.3Y14.013.81

24、4.013.513.512.014.711.413.812.0X, Y分別是治療前,后病人的血紅蛋白的含量數(shù)據(jù),試求治療前后變化的區(qū)間估計.x=c(11.3,15.0,15.0,13.5,12.8,10.0,11.0,12.0,13.0,12.3)y=c(14.0,13.8,14.0,13.5,13.5,12.0,14.7,11.4,13.8,12.0)t.test(x-y)#配對數(shù)據(jù)做差，然后做單樣本t檢驗，其中含有差的 變化的區(qū)間估計 One Sample t-test data: x - y t = -1.3066, df = 9, p-value = 0.2237alternative

25、 hypothesis: true mean is not equal to 0 95 percent confidence interval: -1.8572881 0.4972881 mean of x -0.68 3.方差比的區(qū)間估計方差比的區(qū)間估計1, 2 已知，分別是1, 2的無偏估計，12( ,)F n n22112222/F2/212( ,)1P FFn n 參數(shù) 的置信度為1-的置信區(qū)間2212/222/2121/21221212(/(/,),)Fn nFn n211222()/()/nnnn12211121221222()/()/niiniiXnXn1, 2 未知未知inte

26、rval_var2-function(x,y, mu=c(Inf, Inf), alpha=0.05) n1-length(x); n2-length(y) if (all(muInf) Sx2-1/n1*sum(x-mu1)2); Sy2-1/n2*sum(y-mu2)2) df1-n1; df2-n2 else Sx2-var(x); Sy2-var(y); df1-n1-1; df2-n2-1 r-Sx2/Sy2 a-r/qf(1-alpha/2,df1,df2) b-r/qf(alpha/2,df1,df2) data.frame(rate=r, df1=df1, df2=df2, a

27、=a, b=b)12(1,1)F nn22112222/SFS/2121/22221221221/(1,1)(1,1)SSnFnSSFnn2/212(1,1)1P FFnn 參數(shù) 的置信度為1-的置信區(qū)間2212/222/2121/21221212(/(/,),)Fn nFn n例子例子4.21：已知兩組數(shù)據(jù)：試用兩種方法作方差比的區(qū)間估計.(1)均值已知1, 2 =80.(2)均值未知.a=c(79.98,80.04,80.02,80.04,80.03,80.03,80.04,79.97,80.05,80.03,80.02,80.00,80.02)b=c(80.02,79.94,79.98,

28、79.97,79.97,80.03,79.95,79.97)source(interval_var2.r)interval_var2(a,b,mu=c(80,80) #均值已知1, 2 =80 rate df1 df2 a b0.7326007 13 8 0.1760141 2.482042interval_var2(a,b) rate df1 df2 a b0.5837405 12 7 0.1251097 2.1052694.3.3非正態(tài)總體的區(qū)間估計非正態(tài)總體的區(qū)間估計設(shè)總體X的均值為，方差為 ，X1,X2,Xn為總體X的一個樣本，當(dāng)n充分大時，2interval_estimate3-fu

29、nction(x,sigma=-1,alpha=0.05) n-length(x); xb=0) tmp-sigma/sqrt(n)*qnorm(1-alpha/2) else tmp-sd(x)/sqrt(n)*qnorm(1-alpha/2) data.frame(mean=xb, a=xb-tmp, b=xb+tmp)參數(shù)的置信度為1-的雙側(cè)置信區(qū)間：未知時例例4.21某公司欲估計自己生產(chǎn)的電池壽命，現(xiàn)從其產(chǎn)品中隨機抽取50只電池做壽命試驗(數(shù)據(jù)由計算機產(chǎn)生，服從均值1/r=2.266(單位：100h)的指數(shù)分布).求該公司生產(chǎn)的電池平均壽命的置信度為95%的置信區(qū)間.x=rexp(50

30、,1/2.266)source(interval_estimate3.r)interval_estimate3(x) mean a b1 2.869167 2.255298 3.483036 , 95%的置信區(qū)間4.3.4單側(cè)置信區(qū)間估計單側(cè)置信區(qū)間估計定義4.7:設(shè)X1,X2,Xn是來自總體X的一個樣本, 是包含在總體分布中的未知參數(shù)，對于給定的(0 1),若統(tǒng)計量 滿足則稱隨機區(qū)間 是的置信度為1- 的單側(cè)置信區(qū)間, 為的置信度為1- 的單側(cè)置信下限.),(21_nXXX1),(21_nXXXP),_類似的由單側(cè)置信上限的定義.1/XPZn 參數(shù)的置信度為1-的單側(cè)置信區(qū)間1/XPZn ,

31、ZXn,ZXnnZXnZX(1)1/XPtnSn 參數(shù)的置信度為1-的單側(cè)置信區(qū)間(1)1/XPtnSn (1),ntSnX,(1)tnSXn(1)tXnnS(1)tXnnSR實現(xiàn)：實現(xiàn)：interval_estimate4-function(x, sigma=-1, side=0, alpha=0.05) n-length(x); xb=0) # 已知# side（標(biāo)記），當(dāng)標(biāo)記0時(左側(cè))，按置信上限公式求置信區(qū)間 if (side0) tmp-sigma/sqrt(n)*qnorm(1-alpha) a - -Inf; b 0) tmp-sigma/sqrt(n)*qnorm(1-alpha) a - xb-tmp; b - Inf else tmp - sigma/sqrt(n)*qnorm(1-alpha/2) a - xb-tmp; b