统计建模与R软件第四讲

关于统计建模与R软件第四讲第1页，共38页，2022年，5月20日，19点46分，星期五关于统计量的诱导关系：第2页，共38页，2022年，5月20日，19点46分，星期五两个正态母体诱导的统计量：两个完全不同的正态分布母体诱导F分布具有相同方差的正态分布母体诱导t分布第3页，共38页，2022年，5月20日，19点46分，星期五主要内容4.1矩法4.2极大似然估计4.3估计量的优良性准则4.4区间估计第4页，共38页，2022年，5月20日，19点46分，星期五思想：用样本矩去估计总体矩，总体矩与总体的参数有关，从而得到总体参数的估计。设总体X的分布函数F(x;θ1……θm)中有m个未知参数,假设总体的m阶原点矩存在,n个样本x1……xn

，令总体的k阶原点矩等于样本的k阶原点矩，即4.1矩法……解此方程组得到则称为参数θk的矩法估计量。第5页，共38页，2022年，5月20日，19点46分，星期五一阶，二阶矩法估计参数：更一般的提法为:利用样本的数字特征作为总体的数字特征的估计.例如:无论总体服从什么分布，其均值和方差分别为:解得均值与方差的矩法点估计:第6页，共38页，2022年，5月20日，19点46分，星期五

设总体服从二项分布B(k;p);k,p为未知参数。X1,x2,……,xn是总体X的一个样本，求参数k,p的矩估计。M1是总体均值(一阶原点矩)M2是总体方差(二阶中心矩)解得：第7页，共38页，2022年，5月20日，19点46分，星期五R实现:(1)#N=20,p=0.7,试验次数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#由解析计算给定结果：>N=m1^2/(m1-m2);N#>[1]21.31695>p=(m1-m2)/m1;p#[1]0.6492486第8页，共38页，2022年，5月20日，19点46分，星期五R实现:(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)}第9页，共38页，2022年，5月20日，19点46分，星期五牛顿法：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是列表,返回函数表达式和函数的Jacobi矩阵;x是迭代初值#初始化#把初值记下来#牛顿法:x1=x0-J-1f#index是示性指标，如果等于1表示牛顿法解存在，否则没有解#函数返回一个列表:根,迭代次数,示性指标,函数值第10页，共38页，2022年，5月20日，19点46分，星期五主函数：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第11页，共38页，2022年，5月20日，19点46分，星期五极大似然法定义1：设总体X的概率密度函数或分布律为是未知参数，为来自总体X的样本，称为θ的似然函数(likelihoodfunction).定义2：设总体X的概率密度函数或分布律为是未知参数，为来自总体X的样本，为θ的似然函数,若：是一个统计量，且满足：则称

为θ的极大似然估计.第12页，共38页，2022年，5月20日，19点46分，星期五1.似然函数关于θ连续极值条件，得：似然方程。独立同分布的样本，似然函数具有连乘的形式第13页，共38页，2022年，5月20日，19点46分，星期五例子：正态分布对数似然方程：#multiroot()函数计算#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是似然方程#公式计算>mean(x)[1]0.1273094>sum((x-mean(x))^2)/10[1]1.267102$root[1]0.24807940.9077064第14页，共38页，2022年，5月20日，19点46分，星期五2.似然函数关于θ有间断点设总体X服从区间[a;b]的均匀分布，x=x1;……

;xn为来自总体的一组样本，用极大似然估计法估计参数a;b。似然函数为L(a;b,x)不是a;b的连续函数,其似然方程为:不能求解从极大似然估计的定义出发来求L(a;b,x)的最大值,要使L达到最大，那么b-a应该尽可能的小，但是a不能大于min(x),b不能小于max(x),因此a;b的极大似然估计为:第15页，共38页，2022年，5月20日，19点46分，星期五3.θ是离散参数空间一般地：在鱼塘钓出r条鱼，做上记号，然后再钓出s条，发现有x条有标记第二次钓出的鱼的条数x服从超几何分布：似然函数为L(N;x)=P(X=x)似然函数的比为：将数字带入上式得池塘中鱼的总数为：[500*1000/72]=6944例子：在鱼塘捞出500条鱼，做上记号，然后再捞出1000条，发现有72条有标记，试估计鱼塘所有的鱼有多少？第16页，共38页，2022年，5月20日，19点46分，星期五4.在解似然方差时无法得到解析解，采用数值方法设总体X服从Cauchy分布，其概率密度函数为：其中θ为未知参数.X1,X2,……,Xn是总体X的样本,求θ极大似然估计.Cauchy分布的似然函数为：求导求对数似然方程的解析解是困难的，考虑使用数值方法。1.使用uniroot函数:#参数为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第17页，共38页，2022年，5月20日，19点46分，星期五2.使用optmize()函数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只能求最小，最大问题转化为负的最小问题第18页，共38页，2022年，5月20日，19点46分，星期五关于二项分布的极大似然估计：matlab输出的极大似然估计数值解：x=20.00000.7065fval=210.2846%matlabº¯Êý£ºfunctionf=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主函数：x0=[21,0.5]A=[0,1;0,-1;-1,0]b=[1;0;-20][x,fval]=fmincon(@fg,x0,A,b)矩法估计值：$root[1]20.91589830.6564385$it[1]5第19页，共38页，2022年，5月20日，19点46分，星期五R实现：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输出结果：$par[1]22.03402140.6179089$value[1]209.5277matlab输出的极大似然估计数值解：x=20.00000.7065fval=210.2846结果对比第20页，共38页，2022年，5月20日，19点46分，星期五区间估计：设总体X的分布函数F(x,θ)含有未知参数θ，对于给定值α(0<α<1),若由样本X1,X2,……,Xn确定的两个统计量，和满足:则称随机区间是参数θ的置信度为1-α的置信区间。置信下限置信上限置信度越短越好第21页，共38页，2022年，5月20日，19点46分，星期五3.1一个正态总体的情况3.1.1均值μ的区间估计:

已知时：参数μ的置信度为1-α的双侧置信区间

未知时：参数μ的置信度为1-α的双侧置信区间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)}#默认σ未知#函数返回一个数据框第22页，共38页，2022年，5月20日，19点46分，星期五例子：

4.14某工厂生产的零件长度X被认为服从N(μ,0.04),现从该产品中随机抽取6个，其长度的测量值如下(单位：mm):试求该零件长度的置信系数为0.95的区间估计。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[]置信区间t.test(x):OneSamplet-testdata:xt=162.1555,df=5,p-value=1.692e-10alternativehypothesis:truemeanisnotequalto095percentconfidenceinterval:14.7130015.18700sampleestimates:meanofx14.95假设：H1拒绝H1

(接受H0)的概率几乎所有的统计软件P-value都是这个意思第23页，共38页，2022年，5月20日，19点46分，星期五当μ已知时：3.1.2方差的区间估计参数

的置信度为1-α的双侧置信区间当μ未知时：参数

的置信度为1-α的双侧置信区间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)}#默认μ未知，未知标志=inf#μ已知时，mu<Inf#μ未知时，mu=Inf第24页，共38页，2022年，5月20日，19点46分，星期五例4.16：用区间估计方法估计例4.15的测量误差σ2，分别对均值μ已知(μ=10)和均值未知两种情况进行讨论。####输入数据，调用编好的程序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第25页，共38页，2022年，5月20日，19点46分，星期五4.3.2两个正态总体的情况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的区间估计(置信度为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)}第26页，共38页，2022年，5月20日，19点46分，星期五例子4.17欲比较甲，乙两种棉花品种的优劣，现假设用它们纺出的棉纱强度分别服从N(μ1,2.182)和N(μ2,1.762)，试验者从这两种棉纱中分别抽取样本X1,X2,…,X100和Y1,Y2,…,Y100(数据用计算机模拟，其随机数的均值分别为5.32和5.76),试给出μ1-μ2的置信度为0.95的区间估计。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第27页，共38页，2022年，5月20日，19点46分，星期五例子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)#两个样本方差不同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)第28页，共38页，2022年，5月20日，19点46分，星期五2.配对数据情形下均值差的区间估计编号12345678910X11.315.015.013.512.810.011.012.013.012.3Y14.013.814.013.513.512.014.711.413.812.0X,Y分别是治疗前,后病人的血红蛋白的含量数据,试求治疗前后变化的区间估计.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)#配对数据做差，然后做单样本t检验，其中含有差的变化的区间估计

OneSamplet-testdata:x-yt=-1.3066,df=9,p-value=0.2237alternativehypothesis:truemeanisnotequalto095percentconfidenceinterval:-1.85728810.4972881meanofx-0.68第29页，共38页，2022年，5月20日，19点46分，星期五3.方差比的区间估计μ1,

μ2

已知，分别是σ1,

σ2的无偏估计，参数

的置信度为1-α的置信区间第30页，共38页，2022年，5月20日，19点46分，星期五μ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)}参数

的置信度为1-α的置信区间第31页，共38页，2022年，5月20日，19点46分，星期五例子4.21：已知两组数据：试用两种方法作方差比的区间估计.(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,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

ratedf1df2ab0.73260071380.17601412.482042interval_var2(a,b)

ratedf1df2ab0.58374051270.12510972.105269第32页，共38页，2022年，5月20日，19点46分，星期五4.3.3非正态总体的区间估计设总体X的均值为μ，方差为，X1,X2,…,Xn为总体X的一个样本，当n充分大时，interval_estimate3<-function(x,sigma=-1,alpha=0.05){n<-length(x);xb<-mean(x)if(sigma>=0)tmp<-sigma/sqrt(n)*qnorm(1-alpha/2)elsetmp<-sd(x)/sqrt(n)*qnorm(1-alpha/2)data.frame(mean=xb,a=xb-tmp,b=xb+tmp)}参数μ的置信度为1-α的双侧置信区间：σ未知时第33页，共38页，2022年，5月20日，19点46分，星期五例4.21某公司欲估计自己生产的电池寿命，现从其产品中随机抽取50只电池做寿命试验(数据由计算机产生，服从均值1/r=2.266(单位：100h)的指数分布).求该公司生产的电池平均寿命的置信度为95%的置信区间.x=rexp(50,1/2.266)source("interval_estimate3.r")interval_estimate3(x)

meanab12.8691672.2552983.483036[,]95%的置信区间第34页，共38页，2022年，5月20日，19点46分，星期五4.3.4单侧置信区间估计定义4.7:设X1,X2,…,Xn是来自总体X的一个样本,θ是包含在总体分布中的未知参数，对于给定的α(0<α<1),若统计量满足则称随机区间是θ的置信度为1-α的单侧置信区间,为θ的置信度为1-

α的单侧置信下限.类似的由单侧置信上限的定义.参数μ的置信度为1-α的单侧置信区间参数μ的置信度为