统计建模与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
•此课件下载可自行编辑修改,仅供参考!
感谢您的支持,我们努力做得更好!
谢谢