4-R软件参数估计

4.参数估计安徽师范大学数学计算机科学学院丁新涛关于统计量的诱导关系：22211(),(,)1niiiSXXXNn(0,1)N(0,1)iN22(),()nm2(0,1),()Nn22S2St221()niin22(1)2N/(,)/nFnmm2F()/tnn2Nt2N)1(~1)(ntSnn)1(~222nSnn两个正态母体诱导的统计量：222(,)iN211(,)iN两个完全不同的正态分布母体诱导F分布22(,)iN21(,)iN21112212222211(1)(1,1)(1)nnnSnFnnnSn具有相同方差的正态分布母体诱导t分布122211222122nnnSnSSnn121212()()(2)11tnnSnn主要内容4.1矩法4.2极大似然估计4.3估计量的优良性准则4.4区间估计思想：用样本矩去估计总体矩，总体矩与总体的参数有关，从而得到总体参数的估计。设总体X的分布函数F(x;θ1……θm)中有m个未知参数,假设总体的m阶原点矩存在,n个样本x1……xn，令总体的k阶原点矩等于样本的k阶原点矩，即4.1矩法……解此方程组得到则称为参数θk的矩法估计量。m^1^,...,k^一阶，二阶矩法估计参数：更一般的提法为:利用样本的数字特征作为总体的数字特征的估计.例如:无论总体服从什么分布，其均值和方差分别为:解得均值与方差的矩法点估计:2211()1niiSXXn20(1)01kMkMppp设总体服从二项分布B(k;p);k,p为未知参数。X1,x2,……,xn是总体X的一个样本，求参数k,p的矩估计。ˆˆ,kpM1是总体均值(一阶原点矩)M2是总体方差(二阶中心矩)1Mkp2(1)kppM解得：kR实现:(1)#N=20,p=0.7,试验次数n=100x-rbinom(100,20,0.7);m1=mean(x)m2=sum((x-mean(x))^2)/100m1[1]13.84m2[1]4.8544#由解析计算给定结果：N=m1^2/(m1-m2);N#[1]21.31695p=(m1-m2)/m1;p#[1]0.6492486kR实现:(2)1221(1)fkpMffkppM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)}1122(1)2ffpkkpJffppkkpkp牛顿法：Newtons-function(fun,x,ep=1e-5,it_max=100){index-0;k-1while(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表示牛顿法解存在，否则没有解#函数返回一个列表:根,迭代次数,示性指标,函数值主函数：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.91589830.6564385$it[1]520(1)01kMkMppp极大似然法定义1：设总体X的概率密度函数或分布律为是未知参数，为来自总体X的样本，称(,),fx12,,,nXXX121(;)(;,,,)(,)nniiLxLxxxfxˆ(();)sup(;)LXXLX为θ的似然函数(likelihoodfunction).定义2：设总体X的概率密度函数或分布律为是未知参数，为来自总体X的样本，为θ的似然函数,若：是一个统计量，且满足：(,),fx12,,,nXXX(;)Lx为12ˆˆˆX)=(X,X,,X)n（则称为θ的极大似然估计.ˆ1.似然函数关于θ连续极值条件，得：(;)0iLX似然方程。独立同分布的样本，似然函数具有连乘的形式ln(;)0iLX例子：正态分布对数似然方程：2211ˆ()niiXXn11ˆniiXxn#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.1273094sum((x-mean(x))^2)/10[1]1.267102$root[1]0.24807940.907706422221ln(,;)ln(2)()22inLxx121()0iLFx222241()022iLnFx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服从超几何分布：(;)()(1;)LNxgNLNx()sxxNsrsNCCPXxC似然函数为L(N;x)=P(X=x)似然函数的比为：()()()NsNrNNsrx22()()NrsNrsNrsNxNˆ[]rsNx将数字带入上式得池塘中鱼的总数为：[500*1000/72]=6944例子：在鱼塘捞出500条鱼，做上记号，然后再捞出1000条，发现有72条有标记，试估计鱼塘所有的鱼有多少？4.在解似然方差时无法得到解析解，采用数值方法设总体X服从Cauchy分布，其概率密度函数为：21111()nniix21(;),[1()]fxxx其中θ为未知参数.X1,X2,……,Xn是总体X的样本,求θ极大似然估计.Cauchy分布的似然函数为：1(;)(,)niiLxfx2ln(;)lnln(1())iLxnx求导2101()niiixx求对数似然方程的解析解是困难的，考虑使用数值方法。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-072.使用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))2#ln(;)lnln(1())iLxnxminimum=0.9021objective=254.4463exitflag=1#θ的近似解#-lnL(θ,x)的近似值$minimum[1]1.03418$objective[1]239.4626matlab解#-lnL=min,则lnL=max,#optimize只能求最小，最大问题转化为负的最小问题关于二项分布的极大似然估计：(,;)(1)iiiiXXnXnXLnpXCpp(,;)()(1)maxiiiiXmnXXnXLnpXCpp()(1)iiiXXnXinpxXCpp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]5ln(,;)(ln)iXnLnpXC()ln(1)imnXplniXp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))ln(,;)(ln)iXnLnpXC()ln(1)imnXplniXpR输出结果：$par[1]22.03402140.6179089$value[1]209.5277matlab输出的极大似然估计数值解：x=20.00000.7065fval=210.2846结果对比区间估计：设总体X的分布函数F(x,θ)含有未知参数θ，对于给定值α(0α1),若由样本X1,X2,……,Xn确定的两个统计量，和满足:112ˆ(,,,)nXXX212ˆ(,,,)nXXX112212ˆˆ((,,,)(,,,))1nnPXXXXXX则称随机区间是参数θ的置信度为1-α的置信区间。12ˆˆ(,)置信下限置信上限置信度越短越好3.1一个正态总体的情况3.1.1均值μ的