R与统计推断
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
均匀分布总体U(a,b)中a，b的矩估计分别为：
i 1
a x(1) min( x),
b x( n ) max( x)
P312 2.设总体X~U0,θ)，现从总体中抽取容量为10的 样本，样本值为0.5 ,1.3,0.6,1.7,2.2,1.2,0.8,1.5,2.0,1.6 试求参数θ的极大似然估计值(MLE)。 X (n) 解 θ的极大似然估计为：
x 取总体X~B(1,p),命中率p的矩估计值为 x , 成败率的估计值： 1 x
y<-scan() 1: 1 1 0 1 0 0 1 0 1 1 1 0 1 1 0 1 0 0 1 0 1 0 1 0 0 1 1 0 1 1 0 1 p<-mean(y);r<-p/(1-p);r [1] 1.285714
1 10 p的极大似然估计值 p x / 10, x ak nk . n k 0
a<-c(0:10 );N<-c(0,1,6,7,23,26,21,12,3,1,0) mean<-a*N/100 p.ml_estim<-mean(x)/10;p.ml_estim [1] 0.499
三、双正态总体方差比的置信区间
var.test(x,y)$conf#输出方差比的95%置信区间 单正态总体方差的区间估计与假设检验需自己编写函数。
四、单样本的总体比率p的置信区间
binom.test(x,n,conf=0.95)$conf#输出比率的95%置信区间 或使用包“Hmisc”
Library(Hmisc)#打开包“Hmisc” bincof(x,n,alpha=0.05)$conf#输出比率的95%置信区间
二、两正态总体均值t-检验
t.test(x,y)#一般情况下双侧假设检验 t.test(x,y,ver.equel=TURE)#方差相等时简单公式 先检验：var.test(x,.y)#检验方差齐次，p值>0.05
P351 例6.6.9做区间估计相应的假设检验问题
例3.2.1 一个顾客买了一包标有500g重的一包红糖，觉得
例3.1.4 （P345 例6.6.5） 假设 轮胎寿命服从正态分
布随机抽12只轮胎试用，测得它们的寿命（单位：万千 米）为：4.68 4.85 4.32 4.85 4.62 5.02 5.20 4.60 4.58 4.73 4.38 4.70 .试求平均寿命的0.95的置信区间。 假设轮胎寿命X~N(μ,2) t<-scan(） 4.68 4.85 4.32 4.85 4.62 5.02 5.20 4.60 4.58 4.73 4.38 4.70 t.test(t,conf=0.95)$conf [1] 4.553454 4.868212 #双侧置信区间，置信上、下限 attr(,"conf.level") [1] 0.95 t.test(t,,alternative=“greater”,conf=0.95)$conf [1] 4.58242 Inf#单侧置信区间，置信下限 attr(,"conf.level") [1] 0.95
（三）求（对数）似然函数的极值点， 使用optimize(f,interval,maximun=T),
f<-function(p){logL = 499*log(p)+501*log(1-p)} optimize(f,c(0,1),maximum = TRUE)
$maximum $objective#目标函数在近似解处的函数值 [1] 0.4990094 [1] -693.1452 对于似然方程不易求解的情况，使用（二），（三）.
第三节
3.1 参数估计
R与统计推断
3.1.1 点估计 3.1.2 区间估计 3.1.3 Basyes估计
3.2 假设检验
3.2.1 参数假设检验 3.2.2 正态性检验与非参数检验 3.2.3 分布拟合优度检验 3.2.4 列联表独立性检验
3.1 参数估计
3.1.1点估计
一、矩估计法:样本矩替代总体矩
3.1.2 区间估计
一、单正态总体均值的置信区间 2已知时 2未知时
[ X u1 2 [ X t1 2 S
n , X u1 2 n , X t1 2 S
n] n]
根据置信区间的结果编写R函数 参见[1]《R语言与统计分析》 汤银才 [2]《统计建模与R语言》 薛毅等 R中没有直接的函数做 2已知时置信区间，需自己编写。 利用R中相应假设检验函数给出置信区间 t.test(x,conf=0.95)$conf#输出 2未知时均值95%置信区间
（二）求似然方程的根，使用uniroot(f,interval)
上例中，似然函数与对数似然函数分别为
L( p) C10i p i (1 p)
x x i 1 n 10 xi xi
n
pi1 (1 p)
10 n xi
i 1
n
p 499 (1 p)501
ln L( p) 499ln p 501ln(1 p)
log L( ) 75log(2 ) 18log(1 ) 20log(1 ) 34log
f<-function(r){logL = 75*log(2-r)+18*log(1r)+20*log(1+r)+34*log(r)} > optimize(f,c(0,1),maximum = TRUE) $maximum [1] 0.4815209 $objective [1] 2.518785 >f<-function(r){logL = 75*log(2-r)+18*log(1r)+20*log(1+r)+34*log(r)} > optimize(f,c(0.6,0.7),maximum = TRUE) $maximum [1] 0.6000733 $objective [1] 0.7720295
均匀分布总体U(a,b)中a，b的矩估计分别为：
a X 3Sn ,
b X 3Sn
P312 2.设总体X~U0,θ)，现从总体中抽取容量为10的 样本，样本值为0.5 ,1.3,0.6,1.7,2.2,1.2,0.8,1.5,2.0,1.6 试求参数θ的矩估计值。 由于EX=θ/2,则θ的矩估计量为： 2 X x<-c(0.5 ,1.3,0.6,1.7,2.2,1.2,0.8,1.5,2.0,1.6) theta.point_estim<2*mean(x) ;theta.point_estim [1] 2.68 记录某个篮球运动员在一次比赛中投篮命中与否，观测数据： 11010010111011010010101001101101 试求这个运动员投篮的成败率的矩估计值。
ak (1 , 1 n k ,m ) EX X i k k ( X 1 , n i 1
k
, X n ), k 1, 2,
m.
总体的矩法估计值基本涵盖于总汇表述统计量中. 使用summary(),mean(),var(),sd()等。 例 3.1.1 正态总体N(μ,2)中μ，2的矩估计分别为： 2 1 n n 1 2 2 X , (Xi X ) S n i 1 n
125
例3.1.3 (P317例6.3.7来自百度文库设一次实验可能有4个结果，发生概率
分别为:(2-θ)/4, (1-θ)/4 ,(1+θ)/4 ,θ/4.其中θ在0,1之间。现进行了197 次实验，结果发生的次数分别为：75,18,70,34，求θ的MLE。 解 似然函数与对数似然函数分别为
1 1 18 1 70 34 L( ) ( )75 ( ) ( ) ( ) (2 ) 75 (1 )18 (1 ) 70 34 2 4 4 4 4
二、极大似然估计法：极大似然原理
极大似然估计 使似然函数或对数似然函数达最大值：
L( ) sup L( ) sup f ( xi , )
i 1
n
log L( ) sup log L( ) sup log f ( xi , )
n
（一）总体的极大似然估计值涵盖于总汇描述统计 量中时，使用summary(),mean(),var(),sd()等。 例 3.1.2 正态总体N(μ,2)中μ，2的极大似然估计分别为： 2 1 n n 1 2 2 X , (Xi X ) S n i 1 n
份量不足，于是找到监督部门；于是监督部门就去商店称 了50包红糖（数据在文件夹《从数据到结论》sugar.txt）； 其中均值（平均重量）是498.35g；这的确比500g少，问是 否能够说明厂家生产的这批红糖平均起来不够份量呢？ 于是需要统计检验。可以画出这些重量的直方图 > weight<-scan("sugar.txt"); Read 50 items > hist(weight) 这个直方图看上去象是 正态分布的样本。不妨假 定这一批袋装红糖有正态 分布N(μ,2)。
x<-c(0.5 ,1.3,0.6,1.7,2.2,1.2,0.8,1.5,2.0,1.6) theta.ml_estim<-max(x);theta.ml_estim [1] 2.2
P322 4. 一地质学家为研究密歇根湖的湖滩地区的岩石成分， 取N=100个样品，每个样品中有10块石子，他记录了每个样品 中属石灰石的石子数，假设采样过程相互独立。得到下了数据 0 1 2 3 4 5 6 7 8 9 10 样品中的石子数ak 样品数nk 0 1 6 7 23 26 21 12 3 1 0 求这地区的石子中石灰石的比例p的极大似然估计值(MLE)。 解 设总体X为每个样品中属石灰石的石子数，知：X~B(10,p)
五、两样本的总体比率差的置信区间
prop.test(c(x1,x2),c(n1,n2),conf=0.95)$conf
#输出比率差的95%置信区间
3.2 假设检验
3.2.1 参数假设检验
一、单正态总体均值t-检验
t.test(x,mu=μ0)#双侧假设检验 t.test(x,mu=μ0,alternative=“greater”)#单侧”>” t.test(x,mu=μ0,alt=“less”)#单侧“<”
似然方程：499/p-501/(1-p)=0 f<-function(p){499/p-501/(1-p)} p.ml_estim <- uniroot(f,c(0,1)) p.ml_estim
$root#方程的近似解即估计值 [1] 0.4990002 $f.root#f在近似解处的函数值 [1] -0.0007449476 $iter#迭代次数 [1] 5 $estim.prec#与精确值的误差 [1] 6.103516e-05
对数似然方程为
75 18 20 34 0 2 1 1
> f<-function(r)(-75/(2-r))-18/(1-r)+70/(1+r)+34/r > r.ml_estim <- uniroot(f,c(0,1)) > r.ml_estim $root [1] 0.6067443 $f.root [1] 0.0006437802 $iter [1] 7 $estim.prec [1] 6.103516e-05
二、两个个正态总体均值差的置信区间
一般情形，使用下列语句：
t.test(x,y,conf=0.95)$conf#输出均值差的95%置信区间 若用方差相等的假设，需下列语句： var.test(x,y)$p.value#输出方差齐次检验p-值，p>0.05 t.test(x,y,var=T,conf=0.95)$conf#输出假定方差相等时的结果 若是成对数据，需下列语句： t.test(x,y,paired=T,conf=0.95)$conf#输出成对数据时的结果
i 1
a x(1) min( x),
b x( n ) max( x)
P312 2.设总体X~U0,θ)，现从总体中抽取容量为10的 样本，样本值为0.5 ,1.3,0.6,1.7,2.2,1.2,0.8,1.5,2.0,1.6 试求参数θ的极大似然估计值(MLE)。 X (n) 解 θ的极大似然估计为：
x 取总体X~B(1,p),命中率p的矩估计值为 x , 成败率的估计值： 1 x
y<-scan() 1: 1 1 0 1 0 0 1 0 1 1 1 0 1 1 0 1 0 0 1 0 1 0 1 0 0 1 1 0 1 1 0 1 p<-mean(y);r<-p/(1-p);r [1] 1.285714
1 10 p的极大似然估计值 p x / 10, x ak nk . n k 0
a<-c(0:10 );N<-c(0,1,6,7,23,26,21,12,3,1,0) mean<-a*N/100 p.ml_estim<-mean(x)/10;p.ml_estim [1] 0.499
三、双正态总体方差比的置信区间
var.test(x,y)$conf#输出方差比的95%置信区间 单正态总体方差的区间估计与假设检验需自己编写函数。
四、单样本的总体比率p的置信区间
binom.test(x,n,conf=0.95)$conf#输出比率的95%置信区间 或使用包“Hmisc”
Library(Hmisc)#打开包“Hmisc” bincof(x,n,alpha=0.05)$conf#输出比率的95%置信区间
二、两正态总体均值t-检验
t.test(x,y)#一般情况下双侧假设检验 t.test(x,y,ver.equel=TURE)#方差相等时简单公式 先检验：var.test(x,.y)#检验方差齐次，p值>0.05
P351 例6.6.9做区间估计相应的假设检验问题
例3.2.1 一个顾客买了一包标有500g重的一包红糖，觉得
例3.1.4 （P345 例6.6.5） 假设 轮胎寿命服从正态分
布随机抽12只轮胎试用，测得它们的寿命（单位：万千 米）为：4.68 4.85 4.32 4.85 4.62 5.02 5.20 4.60 4.58 4.73 4.38 4.70 .试求平均寿命的0.95的置信区间。 假设轮胎寿命X~N(μ,2) t<-scan(） 4.68 4.85 4.32 4.85 4.62 5.02 5.20 4.60 4.58 4.73 4.38 4.70 t.test(t,conf=0.95)$conf [1] 4.553454 4.868212 #双侧置信区间，置信上、下限 attr(,"conf.level") [1] 0.95 t.test(t,,alternative=“greater”,conf=0.95)$conf [1] 4.58242 Inf#单侧置信区间，置信下限 attr(,"conf.level") [1] 0.95
（三）求（对数）似然函数的极值点， 使用optimize(f,interval,maximun=T),
f<-function(p){logL = 499*log(p)+501*log(1-p)} optimize(f,c(0,1),maximum = TRUE)
$maximum $objective#目标函数在近似解处的函数值 [1] 0.4990094 [1] -693.1452 对于似然方程不易求解的情况，使用（二），（三）.
第三节
3.1 参数估计
R与统计推断
3.1.1 点估计 3.1.2 区间估计 3.1.3 Basyes估计
3.2 假设检验
3.2.1 参数假设检验 3.2.2 正态性检验与非参数检验 3.2.3 分布拟合优度检验 3.2.4 列联表独立性检验
3.1 参数估计
3.1.1点估计
一、矩估计法:样本矩替代总体矩
3.1.2 区间估计
一、单正态总体均值的置信区间 2已知时 2未知时
[ X u1 2 [ X t1 2 S
n , X u1 2 n , X t1 2 S
n] n]
根据置信区间的结果编写R函数 参见[1]《R语言与统计分析》 汤银才 [2]《统计建模与R语言》 薛毅等 R中没有直接的函数做 2已知时置信区间，需自己编写。 利用R中相应假设检验函数给出置信区间 t.test(x,conf=0.95)$conf#输出 2未知时均值95%置信区间
（二）求似然方程的根，使用uniroot(f,interval)
上例中，似然函数与对数似然函数分别为
L( p) C10i p i (1 p)
x x i 1 n 10 xi xi
n
pi1 (1 p)
10 n xi
i 1
n
p 499 (1 p)501
ln L( p) 499ln p 501ln(1 p)
log L( ) 75log(2 ) 18log(1 ) 20log(1 ) 34log
f<-function(r){logL = 75*log(2-r)+18*log(1r)+20*log(1+r)+34*log(r)} > optimize(f,c(0,1),maximum = TRUE) $maximum [1] 0.4815209 $objective [1] 2.518785 >f<-function(r){logL = 75*log(2-r)+18*log(1r)+20*log(1+r)+34*log(r)} > optimize(f,c(0.6,0.7),maximum = TRUE) $maximum [1] 0.6000733 $objective [1] 0.7720295
均匀分布总体U(a,b)中a，b的矩估计分别为：
a X 3Sn ,
b X 3Sn
P312 2.设总体X~U0,θ)，现从总体中抽取容量为10的 样本，样本值为0.5 ,1.3,0.6,1.7,2.2,1.2,0.8,1.5,2.0,1.6 试求参数θ的矩估计值。 由于EX=θ/2,则θ的矩估计量为： 2 X x<-c(0.5 ,1.3,0.6,1.7,2.2,1.2,0.8,1.5,2.0,1.6) theta.point_estim<2*mean(x) ;theta.point_estim [1] 2.68 记录某个篮球运动员在一次比赛中投篮命中与否，观测数据： 11010010111011010010101001101101 试求这个运动员投篮的成败率的矩估计值。
ak (1 , 1 n k ,m ) EX X i k k ( X 1 , n i 1
k
, X n ), k 1, 2,
m.
总体的矩法估计值基本涵盖于总汇表述统计量中. 使用summary(),mean(),var(),sd()等。 例 3.1.1 正态总体N(μ,2)中μ，2的矩估计分别为： 2 1 n n 1 2 2 X , (Xi X ) S n i 1 n
125
例3.1.3 (P317例6.3.7来自百度文库设一次实验可能有4个结果，发生概率
分别为:(2-θ)/4, (1-θ)/4 ,(1+θ)/4 ,θ/4.其中θ在0,1之间。现进行了197 次实验，结果发生的次数分别为：75,18,70,34，求θ的MLE。 解 似然函数与对数似然函数分别为
1 1 18 1 70 34 L( ) ( )75 ( ) ( ) ( ) (2 ) 75 (1 )18 (1 ) 70 34 2 4 4 4 4
二、极大似然估计法：极大似然原理
极大似然估计 使似然函数或对数似然函数达最大值：
L( ) sup L( ) sup f ( xi , )
i 1
n
log L( ) sup log L( ) sup log f ( xi , )
n
（一）总体的极大似然估计值涵盖于总汇描述统计 量中时，使用summary(),mean(),var(),sd()等。 例 3.1.2 正态总体N(μ,2)中μ，2的极大似然估计分别为： 2 1 n n 1 2 2 X , (Xi X ) S n i 1 n
份量不足，于是找到监督部门；于是监督部门就去商店称 了50包红糖（数据在文件夹《从数据到结论》sugar.txt）； 其中均值（平均重量）是498.35g；这的确比500g少，问是 否能够说明厂家生产的这批红糖平均起来不够份量呢？ 于是需要统计检验。可以画出这些重量的直方图 > weight<-scan("sugar.txt"); Read 50 items > hist(weight) 这个直方图看上去象是 正态分布的样本。不妨假 定这一批袋装红糖有正态 分布N(μ,2)。
x<-c(0.5 ,1.3,0.6,1.7,2.2,1.2,0.8,1.5,2.0,1.6) theta.ml_estim<-max(x);theta.ml_estim [1] 2.2
P322 4. 一地质学家为研究密歇根湖的湖滩地区的岩石成分， 取N=100个样品，每个样品中有10块石子，他记录了每个样品 中属石灰石的石子数，假设采样过程相互独立。得到下了数据 0 1 2 3 4 5 6 7 8 9 10 样品中的石子数ak 样品数nk 0 1 6 7 23 26 21 12 3 1 0 求这地区的石子中石灰石的比例p的极大似然估计值(MLE)。 解 设总体X为每个样品中属石灰石的石子数，知：X~B(10,p)
五、两样本的总体比率差的置信区间
prop.test(c(x1,x2),c(n1,n2),conf=0.95)$conf
#输出比率差的95%置信区间
3.2 假设检验
3.2.1 参数假设检验
一、单正态总体均值t-检验
t.test(x,mu=μ0)#双侧假设检验 t.test(x,mu=μ0,alternative=“greater”)#单侧”>” t.test(x,mu=μ0,alt=“less”)#单侧“<”
似然方程：499/p-501/(1-p)=0 f<-function(p){499/p-501/(1-p)} p.ml_estim <- uniroot(f,c(0,1)) p.ml_estim
$root#方程的近似解即估计值 [1] 0.4990002 $f.root#f在近似解处的函数值 [1] -0.0007449476 $iter#迭代次数 [1] 5 $estim.prec#与精确值的误差 [1] 6.103516e-05
对数似然方程为
75 18 20 34 0 2 1 1
> f<-function(r)(-75/(2-r))-18/(1-r)+70/(1+r)+34/r > r.ml_estim <- uniroot(f,c(0,1)) > r.ml_estim $root [1] 0.6067443 $f.root [1] 0.0006437802 $iter [1] 7 $estim.prec [1] 6.103516e-05
二、两个个正态总体均值差的置信区间
一般情形，使用下列语句：
t.test(x,y,conf=0.95)$conf#输出均值差的95%置信区间 若用方差相等的假设，需下列语句： var.test(x,y)$p.value#输出方差齐次检验p-值，p>0.05 t.test(x,y,var=T,conf=0.95)$conf#输出假定方差相等时的结果 若是成对数据，需下列语句： t.test(x,y,paired=T,conf=0.95)$conf#输出成对数据时的结果