統(tǒng)計(jì)建模與軟件第四講詳解演示文稿

統(tǒng)計(jì)建模與軟件第四講詳解演示文稿當(dāng)前1頁，總共34頁。優(yōu)選統(tǒng)計(jì)建模與軟件第四講當(dāng)前2頁，總共34頁。兩個(gè)正態(tài)母體誘導(dǎo)的統(tǒng)計(jì)量：兩個(gè)完全不同的正態(tài)分布母體誘導(dǎo)F分布具有相同方差的正態(tài)分布母體誘導(dǎo)t分布當(dāng)前3頁，總共34頁。主要內(nèi)容4.1矩法4.2極大似然估計(jì)4.3估計(jì)量的優(yōu)良性準(zhǔn)則4.4區(qū)間估計(jì)當(dāng)前4頁，總共34頁。思想：用樣本矩去估計(jì)總體矩，總體矩與總體的參數(shù)有關(guān)，從而得到總體參數(shù)的估計(jì)。設(shè)總體X的分布函數(shù)F(x;θ1……θm)中有m個(gè)未知參數(shù),假設(shè)總體的m階原點(diǎn)矩存在,n個(gè)樣本x1……xn

，令總體的k階原點(diǎn)矩等于樣本的k階原點(diǎn)矩，即4.1矩法……解此方程組得到則稱為參數(shù)θk的矩法估計(jì)量。當(dāng)前5頁，總共34頁。一階，二階矩法估計(jì)參數(shù)：更一般的提法為:利用樣本的數(shù)字特征作為總體的數(shù)字特征的估計(jì).例如:無論總體服從什么分布，其均值和方差分別為:解得均值與方差的矩法點(diǎn)估計(jì):當(dāng)前6頁，總共34頁。

設(shè)總體服從二項(xiàng)分布B(k;p);k,p為未知參數(shù)。X1,x2,……,xn是總體X的一個(gè)樣本，求參數(shù)k,p的矩估計(jì)。M1是總體均值(一階原點(diǎn)矩)M2是總體方差(二階中心矩)解得：當(dāng)前7頁，總共34頁。R實(shí)現(xiàn):(1)#N=20,p=0.7,試驗(yàn)次數(shù)n=100x<-rbinom(100,20,0.7);m1=mean(x)m2=sum((x-mean(x))^2)/100>m1[1]13.84>m2[1]4.8544#由解析計(jì)算給定結(jié)果：>N=m1^2/(m1-m2);N#>[1]21.31695>p=(m1-m2)/m1;p#[1]0.6492486當(dāng)前8頁，總共34頁。R實(shí)現(xiàn):(2)moment_fun<-function(p){

f<-c(p[1]*p[2]-M1,p[1]*p[2]-p[1]*p[2]^2-M2)J<-matrix(c(p[2],p[1],p[2]-p[2]^2,p[1]-2*p[1]*p[2]),nrow=2,byrow=T)list(f=f,J=J)}當(dāng)前9頁，總共34頁。牛頓法：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<-sqrt((x-x1)%*%(x-x1))if(norm<ep){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ù)返回一個(gè)列表:根,迭代次數(shù),示性指標(biāo),函數(shù)值當(dāng)前10頁，總共34頁。主函數(shù)：x<-rbinom(100,20,0.7);n<-length(x)M1<-mean(x);M2<-(n-1)/n*var(x)source("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(norm<ep){index<-1;break}k<-k+1}obj<-fun(x);list(root=x,it=k,index=index,FunVal=obj$f)}K0,p0$root[1]20.91589830.6564385$it[1]5當(dāng)前11頁，總共34頁。極大似然法定義1：設(shè)總體X的概率密度函數(shù)或分布律為是未知參數(shù)，為來自總體X的樣本，稱為θ的似然函數(shù)(likelihoodfunction).定義2：設(shè)總體X的概率密度函數(shù)或分布律為是未知參數(shù)，為來自總體X的樣本，為θ的似然函數(shù),若：是一個(gè)統(tǒng)計(jì)量，且滿足：則稱

為θ的極大似然估計(jì).當(dāng)前12頁，總共34頁。1.似然函數(shù)關(guān)于θ連續(xù)極值條件，得：似然方程。獨(dú)立同分布的樣本，似然函數(shù)具有連乘的形式當(dāng)前13頁，總共34頁。例子：正態(tài)分布對數(shù)似然方程：#multiroot()函數(shù)計(jì)算#e[1]=\mu,e[2]=\sigma,x=樣本model<-function(e,x){n=length(x)F1=sum(x-e[1]);F2=-n/(e[2])^2+sum((x-[1])^2)/e[2]^4C(F1,F2)}x=rnorm(10)multiroot(f=model,start=c(0,1),x=x)#F1=0,F2=0是似然方程#公式計(jì)算>mean(x)[1]0.1273094>sum((x-mean(x))^2)/10[1]1.267102$root[1]0.24807940.9077064當(dāng)前14頁，總共34頁。2.似然函數(shù)關(guān)于θ有間斷點(diǎn)設(shè)總體X服從區(qū)間[a;b]的均勻分布，x=x1;……

;xn為來自總體的一組樣本，用極大似然估計(jì)法估計(jì)參數(shù)a;b。似然函數(shù)為L(a;b,x)不是a;b的連續(xù)函數(shù),其似然方程為:不能求解從極大似然估計(jì)的定義出發(fā)來求L(a;b,x)的最大值,要使L達(dá)到最大，那么b-a應(yīng)該盡可能的小，但是a不能大于min(x),b不能小于max(x),因此a;b的極大似然估計(jì)為:當(dāng)前15頁，總共34頁。3.θ是離散參數(shù)空間一般地：在魚塘釣出r條魚，做上記號，然后再釣出s條，發(fā)現(xiàn)有x條有標(biāo)記第二次釣出的魚的條數(shù)x服從超幾何分布：似然函數(shù)為L(N;x)=P(X=x)似然函數(shù)的比為：將數(shù)字帶入上式得池塘中魚的總數(shù)為：[500*1000/72]=6944例子：在魚塘撈出500條魚，做上記號，然后再撈出1000條，發(fā)現(xiàn)有72條有標(biāo)記，試估計(jì)魚塘所有的魚有多少？當(dāng)前16頁，總共34頁。4.在解似然方差時(shí)無法得到解析解，采用數(shù)值方法設(shè)總體X服從Cauchy分布，其概率密度函數(shù)為：其中θ為未知參數(shù).X1,X2,……,Xn是總體X的樣本,求θ極大似然估計(jì).Cauchy分布的似然函數(shù)為：求導(dǎo)求對數(shù)似然方程的解析解是困難的，考慮使用數(shù)值方法。1.使用uniroot函數(shù):#參數(shù)為1的cauchy分布x=rcauchy(100,1)f<-function(p)sum((x-p)/(1+(x-p)^2))out<-uniroot(f,c(0,5))>out$root[1]0.9020655$f.root[1]1.800204e-07當(dāng)前17頁，總共34頁。2.使用optmize()函數(shù)x=rcauchy(100,1)loglike<-function(p){n=length(x);log(3.14159)*n+sum(log(1+(x-p)^2))}>optimize(loglike,c(0,5))minimum=0.9021objective=254.4463exitflag=1#θ的近似解#-lnL(θ,x)的近似值$minimum[1]1.03418$objective[1]239.4626matlab解#-lnL=min,則lnL=max,#optimize只能求最小，最大問題轉(zhuǎn)化為負(fù)的最小問題當(dāng)前18頁，總共34頁。關(guān)于二項(xiàng)分布的極大似然估計(jì)：matlab輸出的極大似然估計(jì)數(shù)值解：x=20.00000.7065fval=210.2846%matlaboˉêy￡ofunctionf=fg(sita)x=load('abc.txt');s=0;fori=1:100s1=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.5]A=[0,1;0,-1;-1,0]b=[1;0;-20][x,fval]=fmincon(@fg,x0,A,b)矩法估計(jì)值：$root[1]20.91589830.6564385$it[1]5當(dāng)前19頁，總共34頁。R實(shí)現(xiàn)：obj=function(n){x<-rbinom(100,20,0.7);m<-length(x)f=-sum(log(choose((n[1]%/%1),x)))-(log(n[2])*sum(x)+log(1-n[2])*(m*n[1]-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))R輸出結(jié)果：$par[1]22.03402140.6179089$value[1]209.5277matlab輸出的極大似然估計(jì)數(shù)值解：x=20.00000.7065fval=210.2846結(jié)果對比當(dāng)前20頁，總共34頁。區(qū)間估計(jì)：設(shè)總體X的分布函數(shù)F(x,θ)含有未知參數(shù)θ，對于給定值α(0<α<1),若由樣本X1,X2,……,Xn確定的兩個(gè)統(tǒng)計(jì)量，和滿足:則稱隨機(jī)區(qū)間是參數(shù)θ的置信度為1-α的置信區(qū)間。置信下限置信上限置信度越短越好當(dāng)前21頁，總共34頁。3.1一個(gè)正態(tài)總體的情況均值μ的區(qū)間估計(jì):

已知時(shí)：參數(shù)μ的置信度為1-α的雙側(cè)置信區(qū)間

未知時(shí)：參數(shù)μ的置信度為1-α的雙側(cè)置信區(qū)間interval_estimate1<-function(x,sigma=-1,alpha=0.05){

n<-length(x);xb<-mean(x)if(sigma>=0){tmp<-sigma/sqrt(n)*qnorm(1-alpha/2);df<-n}else{tmp<-sd(x)/sqrt(n)*qt(1-alpha/2,n-1);df<-n-1}data.frame(mean=xb,df=df,a=xb-tmp,b=xb+tmp)}#默認(rèn)σ未知#函數(shù)返回一個(gè)數(shù)據(jù)框當(dāng)前22頁，總共34頁。例子：

4.14某工廠生產(chǎn)的零件長度X被認(rèn)為服從N(μ,0.04),現(xiàn)從該產(chǎn)品中隨機(jī)抽取6個(gè)，其長度的測量值如下(單位：mm):試求該零件長度的置信系數(shù)為0.95的區(qū)間估計(jì)。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)

meandfab14.95614.7899715.11003[]置信區(qū)間t.test(x):OneSamplet-testdata:xt=162.1555,df=5,p-value=1.692e-10alternativehypothesis:truemeanisnotequalto095percentconfidenceinterval:14.7130015.18700sampleestimates:meanofx14.95假設(shè)：H1拒絕H1

(接受H0)的概率幾乎所有的統(tǒng)計(jì)軟件P-value都是這個(gè)意思當(dāng)前23頁，總共34頁。當(dāng)μ已知時(shí)：3.1.2方差的區(qū)間估計(jì)參數(shù)

的置信度為1-α的雙側(cè)置信區(qū)間當(dāng)μ未知時(shí)：參數(shù)

的置信度為1-α的雙側(cè)置信區(qū)間interval_var1<-function(x,mu=Inf,alpha=0.05){

n<-length(x)if(mu<Inf){S2<-sum((x-mu)^2)/n;df<-n}else{S2<-var(x);df<-n-1}a<-df*S2/qchisq(1-alpha/2,df)b<-df*S2/qchisq(alpha/2,df)data.frame(var=S2,df=df,a=a,b=b)}#默認(rèn)μ未知，未知標(biāo)志=inf#μ已知時(shí)，mu<Inf#μ未知時(shí)，mu=Inf當(dāng)前24頁，總共34頁。例4.16：用區(qū)間估計(jì)方法估計(jì)例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)

vardfab0.055100.026851300.1693885interval_var1(x)vardfab0.0583333390.027598510.1944164當(dāng)前25頁，總共34頁。兩個(gè)正態(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<-mean(y)if(all(sigma>=0))#均值差μ1-μ2的區(qū)間估計(jì)(置信度為1-α){tmp<-qnorm(1-alpha/2)*sqrt(sigma[1]^2/n1+sigma[2]^2/n2)df<-n1+n2}else{if(var.equal==TRUE){Sw<-((n1-1)*var(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/(S1^2/n1^2/(n1-1)+S2^2/n2^2/(n2-1))tmp<-qt(1-alpha/2,nu)*sqrt(S1/n1+S2/n2)df<-nu}}data.frame(mean=xb-yb,df=df,a=xb-yb-tmp,b=xb-yb+tmp)}當(dāng)前26頁，總共34頁。例子4.17欲比較甲，乙兩種棉花品種的優(yōu)劣，現(xiàn)假設(shè)用它們紡出的棉紗強(qiáng)度分別服從N(μ1,2.182)和N(μ2,1.762)，試驗(yàn)者從這兩種棉紗中分別抽取樣本X1,X2,…,X100和Y1,Y2,…,Y100(數(shù)據(jù)用計(jì)算機(jī)模擬，其隨機(jī)數(shù)的均值分別為5.32和5.76),試給出μ1-μ2的置信度為0.95的區(qū)間估計(jì)。x=rnorm(100,5.32,2.18)y=rnorm(100,5.76,1.76)source('interval_estimate2.r')interval_estimate2(x,y,sigma=c(2.18,1.76))

meandfab1-0.5325071200-1.0816470.01663265當(dāng)前27頁，總共34頁。例子4.18x=rnorm(12,501.1,2.4)y=rnorm(17,499.7,4.7)source('interval_estimate2.r')interval_estimate2(x,y,var.equal=TRUE)

meandfab0.921158527-1.875633.717947interval_estimate2(x,y)

meandfab0.921158524.46739-1.5942293.436546>t.test(x,y)#兩個(gè)樣本方差不同WelchTwoSamplet-testdata:xandyt=0.7551,df=24.467,p-value=0.4574alternativehypothesis:truedifferenceinmeansisnotequalto095percentconfidenceinterval:-1.5942293.436546sampleestimates:meanofxmeanofy500.8576499.9365t.test(x,y,var.equal=TRUE)當(dāng)前28頁，總共34頁。2.配對數(shù)據(jù)情形下均值差的區(qū)間估計(jì)編號12345678910X11.315.015.013.512.810.011.012.013.012.3Y14.013.814.013.513.512.014.711.413.812.0X,Y分別是治療前,后病人的血紅蛋白的含量數(shù)據(jù),試求治療前后變化的區(qū)間估計(jì).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檢驗(yàn)，其中含有差的變化的區(qū)間估計(jì)

OneSamplet-testdata:x-yt=-1.3066,df=9,p-value=0.2237alternativehypothesis:truemeanisnotequalto095percentconfidenceinterval:-1.85728810.4972881meanofx-0.68當(dāng)前29頁，總共34頁。3.方差比的區(qū)間估計(jì)μ1,

μ2

已知，分別是σ1,

σ2的無偏估計(jì)，參數(shù)

的置信度為1-α的置信區(qū)間當(dāng)前30頁，總共34頁。μ1,

μ2

未知interval_var2<-function(x,y,mu=c(Inf,Inf),alpha=0.05){n1<-length(x);n2<-length(y)if(all(mu<Inf)){Sx2<-1/n1*sum((x-mu[1])^2);Sy2<-1/n2*sum((y-mu[2])^2)df1<-n1;df2<-n2}else{Sx2<-var(x);Sy2<-var(y);df1<-n1-1;df2<-n2-1}r<-Sx2/Sy2a<-r/qf(1-alpha/2,df1,df2)b<-r/qf(alpha/2,df1,df2)data.frame(rate=r,df1=df1,df2=df2,a=a,b=b)}參數(shù)

的置信度為1-α的置信區(qū)間當(dāng)前31頁，總共34頁。例子4.21：已知兩組數(shù)據(jù)：試用兩種方法作方差比的區(qū)間估計(jì).(1)均值已知μ1,

μ2=80.(2)均值未知.a=c(79.98,80.04,80.02,80.04,80.03,80.03,80.0