版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
关于统计建模与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-α的单侧置信区间参数μ的置信度为
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 微控制器项目可行性实施报告
- 2024年现场总线计算机通讯模板项目提案报告模范
- 擦脸毛巾项目可行性实施报告
- 留言板系统课课程设计
- 节水管理宣传标语
- 电子政务机房建设方案
- 空调课程设计选水冷机组
- 棉网市场环境与对策分析
- 烤鱼外卖运营文案策划方案
- 生物制药课程设计
- 读懂中国茶学习通超星期末考试答案章节答案2024年
- DLT5196-2016 火力发电厂石灰石-石膏湿法烟气脱硫系统设计规程
- 新闻采编与制作职业生涯规划
- 部编版九年级上册语文教材知识练习题:字音字形|诗文名句默写|文言文知识|文学常识|传统文化(含答案)
- 《烧(创)伤的急救复苏与麻醉管理》智慧树知到课后章节答案2023年下中国人民解放军总医院第四医学中心
- 贵州大学-学生手册-易班考试答案
- 证据清单模板
- 《合并同类项》优质课一等奖课件
- 农林项目投资实验:财务分析与经济分析综合练习
- 人教版八年级上册英语全册教案(完整版)教学设计含教学反思
- 商业银行IT运维服务项目方案建议书
评论
0/150
提交评论