统计建模与R软件第四讲只是课件.docx
- 文档编号:12399126
- 上传时间:2023-04-18
- 格式:DOCX
- 页数:22
- 大小:1.64MB
统计建模与R软件第四讲只是课件.docx
《统计建模与R软件第四讲只是课件.docx》由会员分享,可在线阅读,更多相关《统计建模与R软件第四讲只是课件.docx(22页珍藏版)》请在冰豆网上搜索。
统计建模与R软件第四讲只是课件
•统计建模与R软件-第四讲・
(2018)
77汁饨俗&)
(日“二座r切+电-2)fl1
孤——
S2十电^^2
具有相同方差的正态分布母体诱导t分布
主要内容
4.1矩法
4.2极大似然估计
4.3估计量的优良性准则
4.4区间估计
4.1矩法
思想:
用样本矩去估计总体矩,总体矩与总体的参数有郑从而得到总体参数的估计。
设总体X的分布函数F(xQ……0m)中有m个未知参数,假设总体的m阶原点矩存注n个样本X[点矩等于样本的k阶原点矩,即
E(X)
=
1n
i=l
1A2二A
i=l
解此方程组得到
✓KZK
则称心为参数乞的矩法估计量。
E(Xm)=
1n
i=l
V
xn,令总体的k阶原
一阶,二阶矩法估计参数
更一般的提法为:
利用样本的数字特征作为总体的数字特征的估计•例如:
无论总体服从什么分布,其均值和方差分别为:
1n
ni=l
E(X)
E(X2)
Var(X)+[E(X)F二/十"
n
i=l
解得均值与方差的矩法点估计:
/X
宀治5
n
丄2—>X,—Xn乙—
i=l
江⑵-疔
i=l
设总体服从二项分布B(k;p);k3p为未知参数。
X15x25……5xn是总体X的一个样本,求参数k,p的矩估计「p;
两—M1=O
\kp(X-p)-M2=Q
Ml2
Ml-M27
M1是总体均值(一阶原点矩)
Ml=kp
M2是总体方差(二阶中心矩)
Ml—M2
R实现:
⑴
#N=205p=0.75试验次数n=100x<-rbinom(100,20,0.7);m1=mean(x)m2=sum((x-mean(x))A2)/100
>ml
[1]13.84
>m2
[1)4.8544
#由解析计算给定结果:
>N=m1A2/(m1-m2);N#1=>[1]21.31695
>p=(m1-m2)/m1;p#p=
[1]0.6492486
Ml2
Ml-M2
Ml—M2
Ml
R实现:
⑵
moment_fun<-function(p){
f<-c(p[irp[2]-M15p[1]*p[2]-p[1]*p[2]A2-M2)
J<-matrix(c(p[2],p[1]5p[2卜p[2]馅p[1]-2*p[1]*p[2])3nrow=23byrow=T)
list(f=f5J=J)
£p」H/Azs@第¥
厂J
-
y
・JcIY#fun是列表,返回函数表达式和
吧exv:
O;k<-1#呷口化函数的Jacobi矩阵;x是迭代初值while(k<=it_max){
x1<-x;#把初值记下来
obj<-fun(x);
x<-x-solve(obj$J5obj$f);#牛顿法:
X[=Xo・J」fnorm<-sqrt((x-x1)%*%(x-x1))
if(norm {index<-1;break}#jndex是示性指标,如果等于1 kv-k+1}表示牛顿法解存在,否则没有解 obj<-fun(x); list(root=x3it=k5index=index5FunVal=obj$f) #函数返回一个列表: 根,迭代次数,示性指标,函数值 obj<-fun(x); 9iist(root二x,it二k,index二index,FunVal=obj$f) } 极大似然法 定义仁设总体X的概率密度函数或分布律为心乡esJ 是未知参数,呂否丄,否为来自总体x的样本,称 为e的似然函数(likelihoodfunction). 定义2: 设总体X的概率密度函数或分布律为心乡名E是未知参数,呂卞丄,氏为来自总体X的样本&朗为 为。 的似然函数,若: 旨是一个统计量,且满足: 则称丨为e的极大似然估计. 1■似然函数关于e连续 极值条件,得: 独立同分布的样本,似然函数具有连乘的形式 F2=-n/(e[2])A2+sum((x-[1])A2)/e[2]A4 C(F1,F2)} x=rnorm(10) multiroot(f=model,start=c(0,1)Jx=x) #F1=0,F2=0是似然方程 $root [1]0.24807940.9077064 2.似然函数关于e有间断点 设总体X服从区间[a;b]的均匀分布,x=x1;……;Xn为》总体的一组样本,用极大似然估计法估计参数a;b。 似然函数为 a<—Xi<—I)、? =1,•…、nothers 从极大似然估计的定义出发来求L(a;b5x)的最大值,要使L达到最大,那么b-a应该尽可能的小,但是a不能大于min(x)3b不能小于max(x),因此a;b的极大似然估计为: d=min(x).b=max(x) 3旧是离散参数空间 例子: 在鱼塘捞出500条鱼,做上记号,然后再捞出10轉条心| 发现有72条有标记,试估计鱼塘所有的鱼有多少? 一般地: 在鱼塘钓出r条鱼,做上记号,然后再钓出s条,发现有x条有标记第二次钓出的鱼的条数x服从超几何分布: p^X=x)= 当rs>2;N时有g(z;N)>1? vs 将数字带入上式得池塘中鱼的总数为: [500*1000/72]=6944 4■在解似然方程时无法得到解析解,采用数值方法 设总体X服从Cauchy分布,其概率密度函数为: 、 其中9为未知参数.X],X2,……,Xn是总体X的样本,求9极大似然估计. Cauchy分布的似然函数为: n1/r1 ■—・-U求导士 求对数似然方程的解析解是困难的, 1•使用uniroot函数: #参数为1的cauchy分布x=rcauchy(10051) f<-function(p)sum((x-p)/(1+(x-p)A2))out<-uniroot(f,c(0,5)) x=rcauchy(1OO,1) loglike<-function(p)#optimize只能求最小,最大问题转化为负的最小问题{n=length(x);_ log(3.14159)*n+sum(log(1+(x-p)A2))/ >optimize(loglike,c(0,5)) $minimum#8的近似解 [1]1.03418 #-lnL=min,贝iJlnL=max, minimum=0.9021matlab解 丄► objective ^objective#-|nL(0,x)的近 =254.4463 exitflag=1 ⑴皿6似值 matlab输出的极大似然估计数值解: x=20.00000.7065 fval=210.2846 R实现: obj=function(n) x<-rbinom(100,20,0.7); m<-length(x)f=・sum(log(choose((n[1]%/%1),x)))二 —|[log(n[5])心um(x)|+|og(1・n[2]广(m*n[1卜sum(x))) sita0=c(20,0.5) #初值 constrOptim(sitaO,obj,NULL,ui=rbind(c(0,- 1),c(O51),c(1,O))3ci=c(-1,0,-20)) R输出结果: $par [1]22.03402140.6179089 $value [1]209.5277 f.屮右11,matlab输出的极大似然估计数值解: ,口果河%=20.00000.7065 fval=210.2846 区间估计: 设总体X的分布函数F(xQ)含有未知参数3对于给定值a(Ovav1),若 本£,X2,人确定的两个统计量,帚化」©和豪応4_,“满足: 3.1.1均值卩的区间估计: 。 : 已知时: y_参数|J的置信度为I-。 的双侧置信区间 —rrE°,l) (J/yjn s/4n interval_estimate1<-function(x|sigma=-l],alpha/zQ.05){n<-length(x);xb<-mean(x)#默认&未知/ if(sigma>=0){else;tmp<-bd(x)/sqrt(n)*qt(1-alpha/2,n-1);df<-n-1data.frame(mean=xb,df=df,a=xb-tmp,b=xb+tmp) }#函数返回一个数据框 例子: 錨6护磋霧龍嚴緇潸隸-冋现从该产品 14.6 15.1 14.9 14.8 15.2 15.1 试求该零件长度的置信系数为0.95的区间估计。 source(1nterval_estimate1・R')x=c(14.6,15.1,14.9,14.8,15.2,15.1)interval_estimate1(x,sigma=0.2) mean 14.95 95percentconfideneeinterval: 14.7130015.18700 sampleestimates: meanofx 14.95 拒绝H1(接受HO) 的概率 假设: H1 3/1.2方差八的区间估计 当|1已知时: ILL.《(爲―“2 斤2■2^Z2 07=1b 参数0'的置信度为1-a的双侧置信区间 当U未知时: 参数。 I的置信度为1-a的双侧置信区间 interval_var1<-function(x,mp^lntai|)ha=0.05){ n<-length(x)#默认》未知,未知标志勻“ if(mu b<-df*S2/qchisq(alpha/2,df) data.frame(var=S2,df=df,a=a,b=b)} 分别对均值U已知(U 用区间估计方法估计例4.15的测量误差a2,均值未知两种情况进行讨论。 ####输入数据,调用编好的程序 •x=c(10.1,10,9.8,10.5,9.7,10.1,9.9,10.2,10.3,9.9) •interval_var1(x5mu=10) vardfab 0.055100.026851300.1693885 interval_var1(x) vardfab 0.0583333390.02759851 0.1944164 interval_estimate2<-function(x,y,sigma=c(-1,-1),var.equal=FALSE,qtlpha=0.05){n1<-length(x);n2<-length(y)xbv・mean(x);yb<-mean(y)if(all(si专ma>=0))#均值釜u厂u2的区间估计(置信度为) {tmp<-qnorm(1-alpha/2)*sqrt(sigma[1]A2/n1+sigma[2]A2/n2)dfv・ else{ if(var.equal==TRUE) Swv*(n1广var(x)+(n2・1广var(y))/(n1+n2・2)tmp<-sqrt(Sw*(1/n1+1/n2))*qt(1-alpha/2,n1+n2-2)dfv・n] else S1<-var(x);S2<-var(y) nu<-fS1/n1+S2/n2)A2/(S1A2/n1A2/(n1-1)+S2A2/n2A2/(n2-1))tmp<-qt(1-alpha/2,nu)*sqrl(S1|/Hl+S2/n2)df<-nu data.frame(mean=xb-yb,df=df,a=xb-yb-tmp,b=xb-yb+tmp)} 欲比较甲,乙两种棉花品种的优劣,现假设用它们纺出的棉纱强度分别 服从N(u「2.180和|\1(u2,1.760,试验者从这两种棉纱中分别抽取样本X2,X100WvY2,…,丫血(数据用计算机模拟,其随机数的均值 分别为5.32和5.76),试给岀口厂u2的置信度为。 ・95的区间估计。 x=rnorm(100,5.32,2.18) y=rnorm(100,5.76,1.76) source。 nterval_estimate2.r,) interval_estimate2(x,y,sigma=c(2.18,1.76)) meandfab 1-0.5325071200-1.0816470.01663265 intervalestimate2(x5y,var.equal=TRUE) intervalestimate2(x,y) i mean df 0.921158524.46739-1.5942293.436546 >t・test(x,y) We.chTWoZZ2^同®(ha: equ: : ;RUE)] data: xandy t=0.7551,df=24.467,p-value=0.4574 alternativehypothesis: truediffereneeinmeansisnotequalto095percentconfideneeinterval: -1.5942293.436546 sampleestimates: meanofxmeanofy 500.8576499.9365 X,Y分别是治疗前,后病人的血红蛋白的含量数据,试求治疗前后变化的区间估计. x=c(11.3,15.0,15.0,13.5,12.8,10.0,11.0,12.0,13.0,12.3) 7=0(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-y •t=-1.3066,df=9,p-value=0.2237 •alternativehypothesis: truemeanisnotequalto0 •95percentconfidenceinterval: -1.85728810.4972881 •meanofx •-0.68 3•方差比的区间估计 m15|J2已知 ni 11In =—x(眾一"1尸能=「y? 、x—“2) 2=12i=l ■分别是5,6的无偏估计, 参数员1云的置信度为1-a的置信区间 a2/2 efj? 話ME) II. y(Xi-“|) 甩%◎'局2%化) M2未知 参数^2/^2的置信度 为1-Q的置信区间 interval_var2<-function(x,y,mu=c(lnf,Inf),alpha=0.05){n1<-length(x);n2<-length(y)if(all(mu Sx2v・1/n广sum((x・mu[1])A2);Sy2<-1/r|2*sum((y-mu[2])A2) df1<-n1;df2<-n24 else{Sx2<-var(x);Sy2<-var(y);d11<-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)} 已知两组数据: A: 79.9880.0480.0280.0480.0380.0380.0479.97 80.0580.0380.0280.0080.02 B: 80.0279.9479.9879.9779.9780.0379.9579.97 试用两种方法作方差比的区间估计•⑴均值已知川,p2=80•⑵均值未知. 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。 nterval_var2.r,) interval_var2(a,b,mu=c(80,80))#均值已知p2=80 ratedf1df2ab 0.73260071380.17601412.482042 interval_var2(a,b) ratedf1df2ab 0.58374051270.12510972.105269 设总体X的均值为p,方差为0: X1,X2,...,Xn为总体X的一个样n充分大时, 土X)—TIRtlN(0,1) 参数|J的置信度为19的双侧置信区间: Q未知时「 亍嗥N/2,X+ ■ X_x+ interval_estimate3<-function(x,sigma=-1,alpha=0.05){ n<-length(x);xbv・mean(x) if(sigma>=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) } 例4・21 某公司欲估计自己生产的电池寿命,现从其产品中随机抽取50只电0寿命试验(数据由计算机产生,服从均值1/r=2.266(单位: 100h)的指2分布)•求该公司生产的电池平均寿命的置信度为95%的置信区间. x=rexp(50,1/2.266)source(,'interval_estimate3.rH)interval_estimate3(x) • >95%的置信区间 meanab •12.869167[2.2552983-483036 4・3・4单侧置信区间估计 定义4.7: 设X1,X2,...5Xn是来自总体X的一个样本,0是包含在总体分的未知参数,对于给定的a(Ovav1),若统计量满足 则称随机区间⑹+①是e的置擔度为a的单侧置信区间,&为e的置信度为a的单侧畫信下限. 类似的由单侧置信上限的定义. 1p参数》的置信度为 1-a的单侧置信区间 =1—G. 乂一金乙严 ^=x- —参数|J的置信度为1),十勺《=xp 1-a的单侧置信区间,l, (F対子1)乂十肩 伽―1) 1) R实现: interval_estimate4<-function(x,sigmaside|=O,alphQ=0X)5){nv・length(x);xb<-mean(x) if(sigma>=0){#o已知 av・xb・tmp;bv・Inf} else{tmpv・sigma/sqrt(n)*qnorm(1・alpha/2) a<-xb-tmp;b<-xb+tmp}#默认side=0,求双侧置信区间df<-n} else{if(side<0){tmp<-sd(x)/sqrt(n)*qt(1-alpha,n-1)a<--Inf;b<-xb+t卩p} elseif(side>0){tmp<-sd(x)/sqrt(n)*qt(1-alpha,n-1) a<-xb-tmp;b<-Inf} else{tmp<-sd(x)/sqrt(n)*qt(1-alpha/2,n-1)口 a<-xb-tmp;b<-xb+tmp}#求双7则置信区间- df<-n-1} data.frame(mean=xb,df=df,a=a,b=b)} 从一批灯泡中随机地抽取5只作寿命试验,测得寿命为: 10501100112012501280 设灯泡寿命服从正态分布,求灯泡寿命平均值的置信度为0.95的单侧置信下限。 x=c(1050,1100,1120,1250,1280) source(Hinterval_estimate4.rH) interval_estimate4(x,side=1) meandfab 1116041064.900Inf #side=O求双侧置信区间 interval_var3<-function(x,mu=lnf,side=0,alpha=0.05){ #side=O默认求双侧置信区间硏=丄丈(兀 nv・length(x) if(mu else{p2v・var(x)~~} if(side<0){a<-0#单侧置信区间: EPb<-df*S2/qchisq(alpha,df)} elseif(side>0){a<-df*S2/qchisq(1-alpha,df)b<-Inf}#a单侧置信区间下限 else{a<-df*S2/qchisq(1-alpha/2,df) b<-df*S2/qchisq(alpha/2,df)} /__1\Q2 data.frame(var=S2,df=df,a=a5b=b) 2_(〃一1)S一汽(—1) 例4・23 求下列10个数据的方差的单侧置信区间上限g=0・05) •x=c(10.1,10,9.8,10.5,9.7,10.1,9.9,10.2,10.3,9.9) •source("interval_var3.rH) •interval_var3(x,side=-1) •vardfab 10.05833333900.1578894 •此课件下载可自行编辑修改,仅供参考! 感谢您的支持,我们努力做得更好! 谢谢
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 统计 建模 软件 第四 只是 课件