4-R软件参数估计.ppt

1、4. 参数估计,安徽师范大学数学计算机科学学院 丁新涛,关于统计量的诱导关系：,两个正态母体诱导的统计量：,两个完全不同的正态分布母体诱导F分布,具有相同方差的正态分布母体诱导t分布,主要内容,4.1 矩法 4.2 极大似然估计 4.3 估计量的优良性准则 4.4 区间估计,思想：用样本矩去估计总体矩，总体矩与总体的参数有关，从而得到总体参数的估计。,设总体X的分布函数F(x;1 m )中有m个未知参数, 假设总体的m阶原点矩存在,n个样本x1 xn ，令总体的k阶原点矩等于样本的k阶原点矩，即,4.1 矩法,解此方程组得到 则称 为参数k的矩法估计量。,一阶，二阶矩法估计参数：,更一般的提法

2、为:利用样本的数字特征作为总体的数字特征的估计.例如:无论总体服从什么分布，其均值和方差分别为:,解得均值与方差的矩法点估计:,设总体服从二项分布B(k; p);k, p为未知参数。X1,x2,xn是总体X的一个样本，求参数k,p的矩估计 。,M1是总体均值(一阶原点矩),M2是总体方差(二阶中心矩),解得：,R实现:(1),# N=20,p=0.7, 试验次数n=100 x-rbinom(100, 20, 0.7); m1=mean(x) m2=sum(x-mean(x)2)/100, m1 1 13.84 m2 1 4.8544,# 由解析计算给定结果： N=m12/(m1-m2); N

3、# 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

4、, 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) M

5、1-mean(x); M2-(n-1)/n*var(x) source(moment_fun.R); source(Newtons.R) p-c(10,0.5); Newtons(moment_fun, p),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 -

6、 fun(x); list(root=x, it=k, index=index, FunVal= obj$f) ,$root 1 20.9158983 0.6564385 $it 1 5,极大似然法,定义1：设总体X的概率密度函数或分布律为 是未知参数， 为来自总体X的样本，称,为的似然函数(likelihood function).,定义2：设总体X的概率密度函数或分布律为 是未知参数， 为来自总体X的样本， 为的似然函数,若： 是一个统计量，且满足：,则称 为的极大似然估计.,1.似然函数关于连续,极值条件，得：,独立同分布的样本，似然函数具有连乘的形式,例子：正态分布,对数似然方程：,#

7、 multiroot()函数计算 # e1=mu, e2=sigma, x=样本 model - function(e,x) n=length(x) F1= sum(x-e1); F2= -n/(e2)2 + sum(x-1)2)/e24 C(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.2480794 0.9077064,2.似然函数关于有间断点,设总体X服从区间a

8、; 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=694

9、4,例子：在鱼塘捞出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-uniroot(f,c(0,5), out $root 1 0.9

10、020655 $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,#-lnL=min,则lnL=max,#optimize只能求最小，最大问题转化为负的最小问题,关于二项分布的极大似然估计：,matlab输出的极大似然估计数值解： x = 20.0000

11、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.

12、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输出的

13、极大似然估计数值解： x = 20.0000 0.7065 fval = 210.2846,区间估计：,设总体X的分布函数F(x,)含有未知参数，对于给定值(0 1),若由样本X1,X2,Xn确定的两个统计量， 和 满足:,则称随机区间 是参数的置信度为1- 的置信区间。,3.1一个正态总体的情况,3.1.1均值的区间估计:,已知时：,未知时：,interval_estimate1=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(mea

14、n=xb, df=df, a=xb-tmp, b=xb+tmp) ,#默认未知,#函数返回一个数据框,例子：,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 mea

15、n is not equal to 0 95 percent confidence interval: 14.71300 15.18700 sample estimates: mean of x 14.95,几乎所有的统计软件P-value都是这个意思,当已知时：,3.1.2方差 的区间估计,当未知时：,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,d

16、f) b-df*S2/qchisq(alpha/2,df) data.frame(var=S2, df=df, a=a, b=b),#已知时， 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

17、 0.05833333 9 0.02759851 0.1944164,4.3.2两个正态总体的情况,interval_estimate2=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)*qt(1-alpha/2,n1+n2-2) df-n1+n2-2 else S1-var(x); S2-var(

18、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-nu data.frame(mean=xb-yb, df=df, a=xb-yb-tmp, b=xb-yb+tmp),例子4.17,欲比较甲，乙两种棉花品种的优劣，现假设用它们纺出的棉纱强度分别服从N(1,2.182)和N(2,1.762)，试验者从这两种棉纱中分别抽取样本X1,X2,X100和Y1,Y2,Y100(数据用计算机模拟，其随机数的均值分别为5.32和5.76),试给出1-2的置信度为0.95的

19、区间估计。,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),mean df a b 1 -0.5325071 200 -1.081647 0.01663265,例子4.18,x=rnorm(12,501.1,2.4) y=rnorm(17,499.7,4.7) source(interval_estimate2.r) interval_estimate2(x,y,var.equal=TRUE),mean df

20、a b 0.9211585 27 -1.87563 3.717947,interval_estimate2(x,y),mean df a b 0.9211585 24.46739 -1.594229 3.436546, t.test(x,y) #两个样本方差不同 Welch Two Sample t-test data: x and y t = 0.7551, df = 24.467, p-value = 0.4574 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence

21、 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.配对数据情形下均值差的区间估计,X, 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) #配对数据做差

22、，然后做单样本t检验，其中含有差的 变化的区间估计,One Sample t-test data: x - y t = -1.3066, df = 9, p-value = 0.2237 alternative 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的无偏估计，,1, 2 未知,interval_var2-function(x,y, mu=c(Inf, Inf),

23、 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) data.frame(rate=r, df1=df1, df2=df2, a=a, b=b),例子4.21：,已知两组数据： 试用两种方法作方差比的区间估

24、计.(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,rate df1 df2 a b 0.7326007 13 8 0.1760141 2.482042,interval_var2(a,b)

25、,rate df1 df2 a b 0.5837405 12 7 0.1251097 2.105269,4.3.3非正态总体的区间估计,设总体X的均值为，方差为 ，X1,X2,Xn为总体X的一个样本，当n充分大时，,interval_estimate3=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)的指数 分布).求该公司生产的电池平均寿命的置信度为95%的置信区间.,x=rexp(50,1/2.266) source(interval_estimate3.r) interval_estimate3(x),mean a b 1 2.869167 2.255298 3.483036,4.3.4单侧置信区间估计,定义4.7:设X1,X2,Xn是来自总体X的一个样本, 是包含在总体分布中的未知参数，对于给定的(0 1),若统计量 满足 则称随机区间 是的置信度为1- 的单侧置