R软件参数估计

1、安徽师范大学数学计算机科学学院安徽师范大学数学计算机科学学院丁新涛丁新涛关于统计量的诱导关系：关于统计量的诱导关系：22211() ,( ,)1niiiSXXXNn (0,1)N(0,1)iN22( ),( )nm2(0,1),( )Nn22S2St221( )niin22(1)2N/( ,)/nF n mm2F( )/t nn2Nt2N) 1(1)(ntSnn) 1(222nSnn两个正态母体诱导的统计量：两个正态母体诱导的统计量：222(,)iN 211(,)iN 两个完全不同的正态分布母体诱导F分布22(,)iN 21(,)iN 211 12212222 211(1)(1,1)(1)nn

2、n SnF nnn Sn具有相同方差的正态分布母体诱导t分布122211222122nnn Sn SSnn121212()()(2)11t nnSnn主要内容主要内容4.1 矩法4.2 极大似然估计4.3 估计量的优良性准则4.4 区间估计思想：用样本矩去估计总体矩，总体矩与总体的参数有关，从而得到总体参数的估计。设总体X的分布函数F(x;1 m )中有m个未知参数, 假设总体的m阶原点矩存在,n个样本x1 xn ，令总体的k阶原点矩等于样本的k阶原点矩，即4.1 4.1 矩法矩法解此方程组得到则称 为参数k的矩法估计量。m1,.,k一阶，二阶矩法估计参数：一阶，二阶矩法估计参数：更一般的提法

3、为:利用样本的数字特征作为总体的数字特征的估计.例如:无论总体服从什么分布，其均值和方差分别为:解得均值与方差的矩法点估计:2211()1niiSXXn20(1)01kMkMppp 设总体服从二项分布B(k; p);k, p为未知参数。X1,x2,xn是总体X的一个样本，求参数k,p的矩估计 。,k pM1是总体均值(一阶原点矩)M2是总体方差(二阶中心矩)1Mkp2(1)kppM 解得：kR实现实现:(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-m

4、2)/m1; p # 1 0.6492486kR实现实现:(2)1221(1)fkpMffkppMmoment_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)1122(1)2ffpkkpJffppkkpkp牛顿法：Newtons-function (fun, x, ep=1e-5, it_max=100) index-0; k-1 while (k=it_max) x1 - x; obj - fun(x);

5、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-len

6、gth(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+

7、1 obj - fun(x); list(root=x, it=k, index=index, FunVal= obj$f)K0,p0$root1 20.9158983 0.6564385$it1 520(1)01kMkMppp极大似然法极大似然法定义1：设总体X的概率密度函数或分布律为是未知参数， 为来自总体X的样本，称( , ),f x 12,nXXX121( ; )( ;,)( , )nniiLxLx xxf x ( ();)sup ( ;)LXXLX为的似然函数(likelihood function).定义2：设总体X的概率密度函数或分布律为是未知参数， 为来自总体X的样本，为的似然

8、函数,若： 是一个统计量，且满足：( , ),f x 12,nXXX( ; )Lx为12X)= (X ,X ,X )n（则称 为的极大似然估计.1.似然函数关于似然函数关于连续连续极值条件，得：( ;)0iLX 似然方程。独立同分布的样本，似然函数具有连乘的形式ln ( ;)0iLX 例子：正态分布例子：正态分布对数似然方程：2211()niiXXn11niiXxn# multiroot()函数计算# e1=mu, e2=sigma, x=样本model mean(x)1 0.1273094 sum(x-mean(x)2)/101 1.267102$root1 0.2480794 0.9077

9、06422221ln ( ,; )ln(2)()22inLxx 121()0iLFx222241()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条鱼，做上记号，然后

10、再钓出s条，发现有x条有标记第二次钓出的鱼的条数x服从超几何分布：(; )()(1; )L N xg NL Nx()s xxN srsNCCP XxC似然函数为L(N;x)=P(X=x)似然函数的比为：()()()Ns NrN Nsrx 22()()Nrs NrsNrs NxNrsNx将数字带入上式得池塘中鱼的总数为：500*1000/72=6944例子：例子：在鱼塘捞出500条鱼，做上记号，然后再捞出1000条，发现有72条有标记，试估计鱼塘所有的鱼有多少？4.在解似然方差时无法得到解析解，采用数值方法在解似然方差时无法得到解析解，采用数值方法设总体X服从Cauchy分布，其概率密度函数为：

11、21111 ()nniix21( ; ),1 () f xxx 其中为未知参数.X1,X2,Xn是总体X的样本,求极大似然估计.Cauchy分布的似然函数为：1( ; )( , )niiLxf x2ln ( ; )lnln(1 () )iLxnx 求导2101 ()niiixx求对数似然方程的解析解是困难的，考虑使用数值方法。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-072.使用optmize()函

12、数x=rcauchy(100,1)loglike optimize(loglike,c(0,5)2# ln ( ; )lnln(1 () )iLxnxminimum =0.9021objective =254.4463exitflag =1#的近似解#-lnL(,x)的近似值$minimum1 1.03418$objective1 239.4626matlab解#-lnL=min,则lnL=max, #optimize只能求最小，最大问题转化为负的最小问题关于二项分布的极大似然估计：关于二项分布的极大似然估计：( , ;)(1)iiiiXXn XnXL n p XCpp ( , ;)()(1)

13、maxiiiiXmnXXnXL n p XCpp ()(1)iiiXXn Xinp xXCppmatlab输出的极大似然估计数值解：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,0

14、b=1;0;-20 x,fval = fmincon(fg,x0,A,b)矩法估计值：$root1 20.9158983 0.6564385$it1 5ln ( , ;)(ln)iXnL n p XC ()ln(1)imnXplniXpR实现：实现：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

15、(0,-1),c(0,1),c(1,0), ci=c(-1,0,-20)ln ( , ;)(ln)iXnL n p XC ()ln(1)imnXplniXpR输出结果：$par1 22.0340214 0.6179089$value1 209.5277matlab输出的极大似然估计数值解：x = 20.0000 0.7065fval = 210.2846结果对比区间估计：区间估计：设总体X的分布函数F(x,)含有未知参数，对于给定值(0 1),若由样本X1,X2,Xn确定的两个统计量， 和 满足:112(,)nXXX212(,)nXXX112212( (,)(,)1nnPXXXXXX 则称随机

16、区间 是参数的置信度为1- 的置信区间。12 ( ,) 置信下限置信上限置信度越短越好3.1一个正态总体的情况一个正态总体的情况3.1.1均值的区间估计:2 已知时：(0,1)/XNn/2/2,XXnZZn/21/XPZn 参数的置信度为1-的双侧置信区间2 未知时：(1)/Xt nSn/2/2(1)(1),SStntnXXnn/21/XPtSn 参数的置信度为1-的双侧置信区间interval_estimate1-function(x,sigma=-1,alpha=0.05) n-length(x); xb=0) tmp-sigma/sqrt(n)*qnorm(1-alpha/2); df-

17、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.115.214.814.915.114.6source(interval_estimate1.R)x=c(14.6,15.1,14.9,14.8,15.2,15.1)interval_est

18、imate1(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)的概率几乎所有的统计软

19、件P-value都是这个意思当已知时：3.1.2方差方差 的区间估计的区间估计22221()( )niiXn22n22/22( )1nPn 参数 的置信度为1-的双侧置信区间22222/21/2,( )( )nnnn当未知时：222(1)(1)nSn222/21/22(1)(1)1),(1)nSnSnn22/22(1)(1)1nSPn 参数 的置信度为1-的双侧置信区间2interval_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); d

20、f - 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 b0.055 10 0.02685

21、130 0.1693885interval_var1(x)var df a b0.05833333 9 0.02759851 0.19441644.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 el

22、se 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

23、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)2.配对数据情形下均值差的区间估计配对数据情形下均值差的区间估计编号12345678910X11.315.015.013

24、.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,

25、 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 3.方差比的区间估计方差比的区间估计1, 2 已知，分别是1, 2的无偏估计，12( ,)F n n22112222/F2/212( ,)1P FFn n 参数 的置信度为1-的置信区间2212/222/2121/21221212(/(/,),)Fn nFn n211222()/()/nnnn1221112

26、1221222()/()/niiniiXnXn1, 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) da

27、ta.frame(rate=r, df1=df1, df2=df2, a=a, b=b)12(1,1)F nn22112222/SFS/2121/22221221221/(1,1)(1,1)SSnFnSSFnn2/212(1,1)1P FFnn 参数 的置信度为1-的置信区间2212/222/2121/21221212(/(/,),)Fn nFn n例子例子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.0

28、2,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.1052694.3.3非正态总体的区间估计非正态总体的区间估计设总体X的均值为，方差为 ，X1,X2,Xn为总体

29、X的一个样本，当n充分大时，2interval_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-的双侧置信区间：未知时例例4.21某公司欲估计自己生产的电池寿命，现从其产品中随机抽取50只电池做寿命试验(数据由计算机产生，服从均值1/r=2.266(单位：100h)的指数分布).

30、求该公司生产的电池平均寿命的置信度为95%的置信区间.x=rexp(50,1/2.266)source(interval_estimate3.r)interval_estimate3(x) mean a b1 2.869167 2.255298 3.483036 , 95%的置信区间4.3.4单侧置信区间估计单侧置信区间估计定义4.7:设X1,X2,Xn是来自总体X的一个样本, 是包含在总体分布中的未知参数，对于给定的(0 1),若统计量 满足则称随机区间 是的置信度为1- 的单侧置信区间, 为的置信度为1- 的单侧置信下限.),(21_nXXX1),(21_nXXXP),_类似的由单侧置信上

31、限的定义.1/XPZn 参数的置信度为1-的单侧置信区间1/XPZn ,ZXn,ZXnnZXnZX(1)1/XPtnSn 参数的置信度为1-的单侧置信区间(1)1/XPtnSn (1),ntSnX,(1)tnSXn(1)tXnnS(1)tXnnSR实现：实现：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