统计建模与R软件第四讲只是课件.docx

1、统计建模与R软件第四讲只是课件统计建模与R软件-第四讲(2018)77汁饨俗&)（日“二座r切+ 电-2） fl 1孤S 2 十电2具有相同方差的正态分布母体诱导t分布主要内容4.1矩法4.2极大似然估计4.3估计量的优良性准则4.4区间估计4.1矩法思想：用样本矩去估计总体矩,总体矩与总体的参数有郑 从而得到总体参数的估计。设总体X的分布函数F(xQ 0m )中有m个未知参数，假设 总体的m阶原点矩存注n个样本X 点矩等于样本的k阶原点矩，即E(X)=1 ni=l1 A 2 二 Ai=l解此方程组得到K ZK则称心为参数乞的矩 法估计量。E(Xm)=1 ni=lVxn ,令总体的k阶原一阶，

2、二阶矩法估计参数更一般的提法为:利用样本的数字特征作为总体的数字特征的 估计例如:无论总体服从什么分布，其均值和方差分别为：1 nn i=lE(X)E(X2)Var(X) + E(X)F 二 /十ni=l解得均值与方差的矩法点估计:/X宀治5n丄2 X, X n乙i=l江-疔i=l设总体服从二项分布B（k; p）;k3 p为未知参数。X15x255xn 是总体X的一个样本，求参数k，p的矩估计p；两M1=Okp（X-p）-M2=QMl2Ml-M27M1是总体均值（一阶原点矩）Ml = kpM2是总体方差（二阶中心矩）Ml M2R实现:# N=205p=0.75 试验次数n=100 x ml1

3、13.84 m21)4.8544#由解析计算给定结果：N=m1A2/(m1-m2); N #1 = 1 21.31695 p=(m1-m2)/m1; p # p =1 0.6492486Ml2Ml - M2Ml M2MlR实现:moment_fun-function(p) f-c(pirp2-M15 p1*p2-p1*p2A2-M2)J-matrix(c(p2, p15 p2 卜 p2馅 p1 -2*p1 *p2)3nrow=23 byrow=T)list(f=f5 J=J)p H / A zs第厂 J-yJ c I Y #fun是列表，返回函数表达式和吧exv：O; k-1 #呷口化 函数的

4、Jacobi矩阵;x是迭代初值 while (k=it_max)x1 -x; #把初值记下来obj - fun(x);x - x - solve(obj$J5 obj$f); #牛顿法:X=XoJf norm - sqrt(x-x1) %*% (x-x1)if (normep) index-1; break #jndex是示性指标，如果等于1kv-k+1 表示牛顿法解存在，否则没有解obj - fun(x);list(root=x3 it=k5 index=index5 FunVal= obj$f)#函数返回一个列表:根，迭代次数, 示性指标,函数值obj - fun(x);9 iist(ro

5、ot二x, it二k, index二index, FunVal= obj$f)极大似然法定义仁 设总体X的概率密度函数或分布律为心乡esJ是未知参数，呂否丄，否为来自总体x的样本，称为e的似然函数(likelihood function).定义2：设总体X的概率密度函数或分布律为心乡名E 是未知参数，呂卞丄，氏为来自总体X的样本&朗为为。的似然函数，若：旨是一个统计量, 且满足：则称丨为e的极大似然估计.1 似然函数关于e连续极值条件，得:独立同分布的样本，似然 函数具有连乘的形式F2= -n/(e2)A2 + sum(x-1)A2)/e2A4C(F1,F2)x=rnorm(10)multir

6、oot(f=model,start=c(0,1)Jx=x)#F1=0,F2=0是似然方程$root1 0.2480794 0.90770642.似然函数关于e有间断点设总体X服从区间a; b的均匀分布，x=x1; ；Xn为 总体的一组样本，用极大似然估计法估计参数a; b。似然函数为a Xi 2；N时有g(z;N) 1? vs zJV时有g(z;N) 1,即似然函 数厶(N; 在N =严附近达到最大即N的极大似然估计为将数字带入上式得池塘中鱼的总数为：500*1000/72=69444在解似然方程时无法得到解析解，采用数值方法设总体X服从Cauchy分布,其概率密度函数为： 、其中9为未知参数

7、.X, X2,Xn是总体X的样本，求9极大似然估计.Cauchy分布的似然函数为：n 1 /r 1 -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 optimize(loglike,c(0,5)$minimum #8的近似解1 1.03418#-lnL=mi n,贝 iJlnL=max, minimum =0.9021 matlab 解丄 object

8、iveobjective #-|nL(0,x)的近=254.4463exitflag =1皿6似值matlab输出的极大似然估计数值解:x = 20.0000 0.7065fval= 210.2846R实现:obj=function(n)x-rbinom(100, 20, 0.7);m-length(x) f=sum(log(choose( n1 %/%1 ),x)二|log(n5)心 um(x)|+|og(1 n2广(m* n1 卜 sum(x)sita0=c(20,0.5)#初值constrOptim(sitaO, obj, NULL, ui=rbind(c(0,-1),c(O51),c(

9、1,O)3ci=c(-1,0,-20)R输出结果：$par1 22.0340214 0.6179089$value1 209.5277,f.屮右11，matlab输出的极大似然估计数值解:，口 果河 % = 20.0000 0.7065fval = 210.2846区间估计:设总体X的分布函数F(xQ)含有未知参数3对于给定值a(Ov av1),若本 ,X2, 人确定的两个统计量，帚化和豪応4_,“满足:3.1.1均值卩的区间估计：。：已知时：y_ 参数|J的置信度为I-。的双侧置信区间rr E，l)(J/yjns/4ninterval_estimate1-function(x|sigma=-

10、l,alpha/zQ.05) n -length(x);xb=0) else; tmp-bd(x)/sqrt(n)*qt(1-alpha/2,n-1); df-n-1 data.frame(mean=xb, df=df, a=xb-tmp, b=xb+tmp) #函数返回一个数据框例子:錨6护磋霧龍嚴緇潸隸-冋现从该产品14.615.114.914.815.215.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

11、=0.2)mean14. 9595 percent con fide nee interval:14.71300 15.18700sample estimates:mean of x14.95拒绝H1 （接受HO）的概率假设：H13/1.2方差八的区间估计当|1已知时：ILL.(爲“2斤 2 2 Z20 7=1 b参数0 的置信度为1-a 的双侧置信区间当U未知时:参数。I的置信度为1-a 的双侧置信区间interval_var1-function(x,mplntai|)ha=0.05)n-length(x) #默认未知,未知标志勻“if (mulnf) S2 - sum(x-mu)A2)/n

12、; df v n #|J已知时，muvlnf else S2 - var(x); df - n-1 #|J未知时，mu=lnf 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)分别对均值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)in terval_var1 (x5mu=10)var df a b0.055

13、 10 0.02685130 0.1693885interval_var1 (x)var df a b0.05833333 9 0.027598510.1944164interval_estimate2-function(x, y, sigma=c(-1,-1 ),var.equal=FALSE, qtlpha=0.05) n 1-length(x); n2-length(y) xbvmean(x); yb=0) #均值釜u厂 u 2的区间估计(置信度为) tmp -qnorm(1 -alpha/2)*sqrt(sigma1 A2/n1 +sigma2A2/n2) dfvelseif (var

14、.equal = TRUE)Swv*(n1 广 var(x)+(n21 广 var(y)/(n1 +n 22) tmp-sqrt(Sw*(1/n1 +1/n2)*qt(1 -alpha/2,n1+n2-2) dfvn elseS1-var(x); S2-var(y)nu-fS1/n1+S2/n2)A2/(S1 A2/n1 A2/(n1-1 )+S2A2/n2A2/(n2-1) tmp-qt(1 -alpha/2, nu)*sqrl(S1|/Hl +S2/n2) df ttest(x,y)We.chTWoZZ2同 (ha：equ：；RUE)data: x and yt = 0.7551, df

15、= 24.467, p-value = 0.4574alter native hypothesis: true differe nee in means is not equal to 0 95 perce nt con fide nee interval:-1.594229 3.436546sample estimates: mean of x mean of y500.8576 499.9365X,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.

16、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 - yt = -1.3066, df = 9, p-value = 0.2237alter native hypothesis: true mean is not equal to 095 percent confidence interval: -1.8572881 0.4972881mean of x-0.683 方差比的区间估计m15 |J2 已知ni1 1 I

17、 n=x(眾 一 1尸能=y?、x“2)2=1 2 i=l分别是5, 6的无偏估计，参数员1云的置信度 为1-a的置信区间a 2 / 2efj?話 ME)II .y( X i - “ |)甩局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(mulnf) Sx2v1/n广sum(xmu1 )A2); Sy2-1 /r|2*sum(y-mu2)A2)df1-n1; df2-n2 4else Sx2-var(x);

18、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.98 80.04 80.02 80.04 80.03 80.03 80.04 79.9780.05 80.03 80.02 80.00 80.02B:80.02 79.94 79.98 79.97 79.97 80.03 79.95 79.97试用两种方法作方差比的区间估计均值已知川，p2 =80均

19、值未知.a=c(7998,80.04,80.02,80.04,80.03,80.03,80.04,79.97,80.05,80.03,8 0.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 =80rate df1 df2 a b0.7326007 13 8 0.1760141 2.482042interval_var2(a,b)rate df1 df2 a b0.5837405 12 7

20、0.1251097 2.105269设总体X的均值为p,方差为0 : , X1,X2,.,Xn为总体X的一个样 n充分大时，土 X)TIR t l N(0,1)参数|J的置信度为1 9的双侧置信区间： Q未知时亍嗥N/2, X +X _ x +interval_estimate3-function(x,sigma=-1,alpha=0.05)n=0)tmp-sigma/sqrt(n)*qnorm(1-alpha/2)elsetmp 95%的置信区间mea n a b1 2.869167 2.255298 3-483036434单侧置信区间估计定义4.7:设X1,X2,.5Xn是来自总体X的一个

21、样本,0是包含在总体分 的未知参数，对于给定的a（Ov av1）,若统计量 满足则称随机区间+是e的置擔度为a的单侧置信区间，&为e的置信度 为a的单侧畫信下限.类似的由单侧置信上限的定义.1p参数的置信度为1-a的单侧置信区间=1G.乂一金乙严=x-参数|J的置信度为1），十勺=xp1-a的单侧置信区间 ，l ,（F対子1） 乂十肩伽1)1)R实现:interval_estimate4-function(x, sigmaside|=O, alphQ=0X)5) nvlength(x); xb=0) # o已知a v xbtmp; b v Inf else tmp v sigma/sqrt(n

22、)*qnorm(1 alpha/2)a - xb-tmp; b - xb+tmp #默认side=0,求双侧置信区间 df-n else if (side0) tmp - sd(x)/sqrt(n)*qt(1 -alpha,n-1) a - -Inf; b 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,

23、 a=a, b=b)从一批灯泡中随机地抽取5只作寿命试验，测得寿命为:1050 1100 1120 1250 1280设灯泡寿命服从正态分布，求灯泡寿命平均值的置信度为0.95的单侧置 信下限。x=c(1050,1100,1120,1250,1280)source(Hinterval_estimate4.rH)interval_estimate4(x,side=1)mean df a b1 1160 4 1064.900 Inf#side=O求双侧置信区间interval_var3-function(x,mu=lnf,side=0,alpha=0.05)#side=O默认求双侧置信区间 硏=丄

24、丈(兀nvlen gth(x)if (mulnf) else p2vvar(x) if (side0) a -0 #单侧置信区间:EP b 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 Q 2data.frame(var=S2, df=df, a=a5 b=b)2 _ (一 1)S 一汽(1)例423求下列10个数据的方差的单侧置信区间上限g=005)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) var df a b1 0.05833333 9 0 0.1578894此课件下载可自行编辑修改，仅供参考！ 感谢您的支持，我们努力做得更好！谢谢