《统计建模与r软件第四讲.ppt》由会员分享，可在线阅读，更多相关《统计建模与r软件第四讲.ppt（38页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、关于统计建模与关于统计建模与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，令
2、总体的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
3、,试验次数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
4、)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
5、页，共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
6、 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.似然函数关于似然函数关于连续连续极值条件，得：似然
7、方程。独立同分布的样本，似然函数具有连乘的形式现在学习的是第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的连续函数,其似然方程为:不能求解从
8、极大似然估计的定义出发来求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页，共
9、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=r
10、cauchy(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%matlabfunctio
11、n 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-rbi
12、nom(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页区
13、间估计：区间估计：设总体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);
14、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
15、)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)的概率几乎所有的统计软件
16、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)#默认未知，未知标志
17、=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
18、.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)
19、*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
20、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分别是治疗前,后病人的血红蛋白的含量数据
21、,试求治疗前后变化的区间估计.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 inte
22、rval:-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)
23、;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)
24、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为
25、总体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(单位
26、：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-的单侧置信下限.类似的由单侧置
27、信上限的定义.参数的置信度为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(
28、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页