统计建模与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（mub<-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（muSx2v・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（muelse{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
•此课件下载可自行编辑修改，仅供参考！
感谢您的支持，我们努力做得更好！
谢谢
升级会员 