R软件参数估计资料

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,4.参数估量,安徽师范大学数学计算机科学学院,丁新涛,关于统计量的诱导关系：,两个正态母体诱导的统计量：,两个完全不同的正态分布母体诱导,F,分布,具有一样方差的正态分布母体诱导t分布,主要内容,4.1 矩法,4.2 极大似然估量,4.3 估量量的优良性准则,4.4 区间估量,思想：用样本矩去估量总体矩，总体矩与总体的参数有关，从而得到总体参数的估量。,设总体,X,的分布函数,F(x;,1,m,),中有,m,个未知参数,假设总体的,m,阶原点矩存在,n,个样本,x,1,x,n,，令总体的,k,阶原点矩等于样本的,k,阶原点矩，即,4.1,矩法,解此方程组得到,则称 为参数,k,的矩法估计量。,一阶，二阶矩法估量参数：,更一般的提法为:利用样本的数字特征作为总体的数字特征的估量.例如:无论总体听从什么分布，其均值和方差分别为:,解得均值与方差的矩法点估量:,设总体服从二项分布B,(k,;p,);k,p,为未知参数,。,X1,x2,xn,是总体,X,的一个样本，求参数,k,p,的矩估计,。,M1,是总体均值,(,一阶原点矩,),M2,是总体方差,(,二阶中心矩,),解得：,R,实现,:(1),#N=20,p=0.7,试验次数,n=100,x m1,1 13.84,m2,1 4.8544,#,由解析计算给定结果：,N=m12/(m1-m2);N#,1 21.31695,p=(m1-m2)/m1;p#,1 0.6492486,R,实现,:(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),牛顿法：,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,是迭代初值,#,初始化,#把初值登记来,#,牛顿法,:x,1,=x,0,-J,-1,f,#index是示性指标，假设等于1,表示牛顿法解存在，否则没有解,#,函数返回一个列表,:,根,迭代次数,示性指标,函数值,主函数：,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,J,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),K0,p0,$root,1 20.9158983 0.6564385,$it,1 5,极大似然法,定义,1,：设总体,X,的概率密度函数或分布律为,是未知参数，为来自总体,X,的样本，称,为,的似然函数,(likelihood function).,定义,2,：设总体,X,的概率密度函数或分布律为,是未知参数，为来自总体,X,的样本，,为,的似然函数,若：是一个统计量，且满足：,则称,为,的极大似然估计,.,1.,似然函数关于,连续,极值条件，得：,似然方程。,独立同分布的样本，似然函数具有连乘的形式,例子：正态分布,对数似然方程：,#multiroot函数计算,#e1=mu,e2=sigma,x=样本,model mean(x),1 0.1273094,sum(x-mean(x)2)/10,1 1.267102,$root,1 0.2480794 0.9077064,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的极大似然估量为:,3.,是离散参数空间,一般地：在鱼塘钓出r条鱼，做上记号，然后再钓出s条，觉察有x条有标记其次次钓出的鱼的条数x听从超几何分布：,似然函数为,L(N;x)=P(X=x),似然函数的比为：,将数字带入上式得池塘中鱼的总数为：,500*1000/72=6944,例子：在鱼塘捞出500条鱼，做上记号，然后再捞出1000条，,觉察有72条有标记，试估量鱼塘全部的鱼有多少？,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,$root,1 0.9020655,$f.root,1 1.800204e-07,2.使用optmize函数,x=rcauchy(100,1),loglike optimize(loglike,c(0,5),minimum,=0.9021,objective,=254.4463,exitflag=1,#,的近似解,#-,lnL(,x),的近,似值,$minimum,1,1.03418,$objective,1,239.4626,matlab,解,#-lnL=min,则,lnL=max,#optimize,只能求最小，最大问题转化为负的最小问题,关于二项分布的极大似然估量：,matlab输出的极大似然估量数值解：,x=20.0000 0.7065,fval=210.2846,%matlab,function 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;,end,f=-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.9158983 0.6564385,$it,1 5,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,输出结果：,$par,1 22.0340214 0.6179089,$value,1 209.5277,matlab输出的极大似然估量数值解：,x=20.0000 0.7065,fval=210.2846,结果对比,区间估量：,设总体,X,的分布函数,F(x,),含有未知参数,，对于给定值,(0,1),若由样本,X,1,X,2,X,n,确定的两个统计量，和 满足,:,则称随机区间 是参数,的置信度为,1-,的置信区间。,置信下限,置信上限,置信度,越短越好,3.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),#,默认,未知,#,函数返回一个数据框,例子：,4.14,某工厂生产的零件长度,X,被认为服从,N(,0.04),现从该产品中随机抽取,6,个，其长度的测量值如下,(,单位：,mm):,试求该零件长度的置信系数为,0.95,的区间估计。,15.1,15.2,14.8,14.9,15.1,14.6,source(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-test,data:x,t=162.1555,df=5,p-value=1.692e-10,alternative 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都是这个意思,当时：,3.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,例,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),var df a b,0.055 10 0.02685130 0.1693885,interval_var1(x),var df a b,0.