统计建模与R软件第四讲精选PPT.ppt

关于统计建模与关于统计建模与R R软软件第四讲件第四讲第1页，讲稿共38张，创作于星期二关于统计量的诱导关系：关于统计量的诱导关系：第2页，讲稿共38张，创作于星期二两个正态母体诱导的统计量：两个正态母体诱导的统计量：两个完全不同的正态分布母体诱导F分布具有相同方差的正态分布母体诱导t分布第3页，讲稿共38张，创作于星期二主要内容主要内容4.1 矩法4.2 极大似然估计4.3 估计量的优良性准则4.4 区间估计第4页，讲稿共38张，创作于星期二思想：用样本矩去估计总体矩，总体矩与总体的参数有关，从而得到总体参数的估计。设总体X的分布函数F(x;1 m)中有m个未知参数,假设总体的m阶原点矩存在,n个样本x1 xn，令总体的k阶原点矩等于样本的k阶原点矩，即4.1 4.1 矩法矩法解此方程组得到则称 为参数k的矩法估计量。第5页，讲稿共38张，创作于星期二一阶，二阶矩法估计参数：一阶，二阶矩法估计参数：更一般的提法为:利用样本的数字特征作为总体的数字特征的估计.例如:无论总体服从什么分布，其均值和方差分别为:解得均值与方差的矩法点估计:第6页，讲稿共38张，创作于星期二 设总体服从二项分布B(k;p);k,p为未知参数。X1,x2,xn是总体X的一个样本，求参数k,p的矩估计 。M1是总体均值(一阶原点矩)M2是总体方差(二阶中心矩)解得：第7页，讲稿共38张，创作于星期二R R实现实现:(1):(1)#N=20,p=0.7,试验次数n=100 x m11 13.84 m21 4.8544#由解析计算给定结果：N=m12/(m1-m2);N#1 21.31695 p=(m1-m2)/m1;p#1 0.6492486第8页，讲稿共38张，创作于星期二R R实现实现:(2):(2)moment_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)第9页，讲稿共38张，创作于星期二牛顿法：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(normep)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张，创作于星期二主函数：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(normep)index-1;break k-k+1 obj-fun(x);list(root=x,it=k,index=index,FunVal=obj$f)K0,p0$root1 20.9158983 0.6564385$it1 5第11页，讲稿共38张，创作于星期二极大似然法极大似然法定义1：设总体X的概率密度函数或分布律为是未知参数，为来自总体X的样本，称为的似然函数(likelihood function).定义2：设总体X的概率密度函数或分布律为是未知参数，为来自总体X的样本，为的似然函数,若：是一个统计量，且满足：则称 为的极大似然估计.第12页，讲稿共38张，创作于星期二1.1.似然函数关于似然函数关于连续连续极值条件，得：似然方程。独立同分布的样本，似然函数具有连乘的形式第13页，讲稿共38张，创作于星期二例子：正态分布例子：正态分布对数似然方程：#multiroot()函数计算#e1=mu,e2=sigma,x=样本model mean(x)1 0.1273094 sum(x-mean(x)2)/101 1.267102$root1 0.2480794 0.9077064第14页，讲稿共38张，创作于星期二2.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张，创作于星期二3.3.是离散参数空间是离散参数空间一般地：在鱼塘钓出r条鱼，做上记号，然后再钓出s条，发现有x条有标记第二次钓出的鱼的条数x服从超几何分布：似然函数为L(N;x)=P(X=x)似然函数的比为：将数字带入上式得池塘中鱼的总数为：500*1000/72=6944例子：例子：在鱼塘捞出500条鱼，做上记号，然后再捞出1000条，发现有72条有标记，试估计鱼塘所有的鱼有多少？第16页，讲稿共38张，创作于星期二4.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 out$root1 0.9020655$f.root1 1.800204e-07第17页，讲稿共38张，创作于星期二2.使用optmize()函数x=rcauchy(100,1)loglike optimize(loglike,c(0,5)minimum=0.9021objective=254.4463exitflag=1#的近似解#-lnL(,x)的近似值$minimum1 1.03418$objective1 239.4626matlab解#-lnL=min,则lnL=max,#optimize只能求最小，最大问题转化为负的最小问题第18页，讲稿共38张，创作于星期二关于二项分布的极大似然估计：关于二项分布的极大似然估计：matlab输出的极大似然估计数值解：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主函数：x0=21,0.5A=0,1;0,-1;-1,0b=1;0;-20 x,fval=fmincon(fg,x0,A,b)矩法估计值：$root1 20.9158983 0.6564385$it1 5第19页，讲稿共38张，创作于星期二R实现：实现：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)R输出结果：$par1 22.0340214 0.6179089$value1 209.5277matlab输出的极大似然估计数值解：x=20.0000 0.7065fval=210.2846结果对比第20页，讲稿共38张，创作于星期二区间估计：区间估计：设总体X的分布函数F(x,)含有未知参数，对于给定值(0 1),若由样本X1,X2,Xn确定的两个统计量，和 满足:则称随机区间 是参数的置信度为1-的置信区间。置信下限置信上限置信度越短越好第21页，讲稿共38张，创作于星期二3.13.1一个正态总体的情况一个正态总体的情况3.1.1均值的区间估计:已知时：参数的置信度为1-的双侧置信区间 未知时：参数的置信度为1-的双侧置信区间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/2,n-1);df-n-1 data.frame(mean=xb,df=df,a=xb-tmp,b=xb+tmp)#默认未知#函数返回一个数据框第22页，讲稿共38张，创作于星期二例子：例子：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)mean df a b 14.95 6 14.78997 15.11003 置信区间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假设：H1拒绝H1(接受H0)的概率几乎所有的统计软件P-value都是这个意思第23页，讲稿共38张，创作于星期二当已知时：3.1.23.1.2方差方差 的区间估计的区间估计参数 的置信度为1-的双侧置信区间当未知时：参数 的置信度为1-的双侧置信区间interval_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)b-df*S2/qchisq(alpha/2,df)data.frame(var=S2,df=df,a=a,b=b)#默认未知，未知标志=inf#已知时，muInf#未知时，mu=Inf第24页，讲稿共38张，创作于星期二例例4.164.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)var df a b0.055 10 0.02685130 0.1693885interval_var1(x)var df a b0.05833333 9 0.02759851 0.1944164第25页，讲稿共38张，创作于星期二4.3.24.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=0)#均值差1-2的区间估计(置信度为1-)tmp-qnorm(1-alpha/2)*sqrt(sigma12/n1+sigma22/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/(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.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)第28页，讲稿共38张，创作于星期二2.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检验，其中含有差的 变化的区间估计 One Sample t-test data:x-y t=-1.3066,df=9,p-value=0.2237alternative hypothesis:true mean is not equal to 0 95 percent confidence interval:-1.8572881 0.4972881 mean of x -0.68 第29页，讲稿共38张，创作于星期二3.3.方差比的区间估计方差比的区间估计1,2 已知，分别是1,2的无偏估计，参数 的置信度为1-的置信区间第30页，讲稿共38张，创作于星期二1 1,2 2 未知未知interval_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=a,b=b)参数 的置信度为1-的置信区间第31页，讲稿共38张，创作于星期二例子例子4.214.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 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.105269第32页，讲稿共38张，创作于星期二4.3.34.3.3非正态总体的区间估计非正态总体的区间估计设总体X的均值为，方差为 ，X1,X2,Xn为总体X的一个样本，当n充分大时，interval_estimate3-function(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)参数的置信度为1-的双侧置信区间：未知时第33页，讲稿共38张，创作于星期二例例4.214.21某公司欲估计自己生产的电池寿命，现从其产品中随机抽取50只电池做寿命试验(数据由计算机产生，服从均值1/r=2.266(单位：100h)的指数分布).求该公司生产的电池平均寿命的置信度为95%的置信区间.x=rexp(50,1/2.266)source(interval_estimate3.r)interval_estimate3(x)mean a b1 2.869167 2.255298 3.483036 ,95%的置信区间第34页，讲稿共38张，创作于星期二4.3.44.3.4单侧置信区间估计单侧置信区间估计定义4.7:设X1,X2,Xn是来自总体X的一个样本,是包含在总体分布中的未知参数，对于给定的(0 1),若统计量 满足则称随机区间 是的置信度为1-的单侧置信区间,为的置信度为1-的单侧置信下限.类似的由单侧置信上限的定义.参数的置信度为1-的单侧置信区间参数的置信度为1-的单侧置信区间第35页，讲稿共38张，创作于星期二R R实现：实现：interval_estimate4-function(x,sigma=-1,side=0,alpha=0.05)n-length(x);xb=0)#已知#side（标记），当标记0时(左侧)，按置信上限公式求置信区间 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-xb+tmp#默认side=0,求双侧置信区间 df-n else if(side0)tmp-sd(x)/sqrt(n)*qt(1-alpha,n-1)a-Inf;b 0)tmp-sd(x)/sqrt(n)*qt(1-alpha,n-1)a-xb-tmp;b-Inf else tmp-sd(x)/sqrt(n)*qt(1-alpha/2,n-1)a-xb-tmp;b-xb+tmp#求双侧置信区间 df-n-1 data.frame(mean=xb,df=df,a=a,b=b)第36页，讲稿共38张，创作于星期二例例4.224.22从一批灯泡中随机地抽取5只作寿命试验，测得寿命为：1050 1100 1120 1250 1280设灯泡寿命服从正态分布，求灯泡寿命平均值的置信度为0.95的单侧置信下限。x=c(1050,1100,1120,1250,1280)source(interval_estimate4.r)interval_estimate4(x,side=1)mean df a b1 1160 4 1064.900 Inf第37页，讲稿共38张，创作于星期二感感谢谢大大家家观观看看第38页，讲稿共38张，创作于星期二