第二章 统计推断
2.1参数估计
点估计是用样本统计量 θ ^ \hat{\theta} θ^的某个取值直接作为总体参数的估计值，因为统计量包含了数据信息，具有数据特征，比如样本均值或者样本比例，也称为矩估计。但是在用矩估计法时，具体样本还需要围绕矩估计构建区间估计，用于概率度量真实值与估计值之间的接近程度。
2.1.2 单总体参数区间估计
##3σ
pnorm(1)-pnorm(-1)
pnorm(2)-pnorm(-2)
pnorm(3)-pnorm(-3)
两种说法：大约有95%的样本均值会落在μ的两个标准差以内，或说，约有95%的样本所构造的两个标准差区间会包括μ。（90%临界值1.645；95%临界值1.96；99%置信水平的临界值2.58）
#总体比例的区间估计
##基于索尔克的试验
n<-200745
(p.hat<-57/n)
p.hat+c(-1.96,1.96)*sqrt(p.hat*(1-p.hat)/n)
##Wald方法是利用正态分布来对二项分布近似，Clopper-Person是完全基于二项分布，是更确切的区间估计方法
binom.test(57,200745)
总体均值的区间估计
##方差已知
#写一个函数用于，方差已知时置信区间的计算函数
conf.inf<-function(x,n,sigma,alpha){
options(digits=5)
mean<-mean(x)
c(mean-sigma*qnorm(1-alpha/2,mean=0,sd=1,
lower.tail=TRUE)/sqrt(n),
mean+sigma*qnorm(1-alpha/2,mean=0,sd=1,
lower.tail=TRUE)/sqrt(n))
}
x<-c(112.5,101.0,...,93.3)#数据集读入
n=25;alpha=0.05;sigma=10;
result<-conf.int(x,n,sigma,alpha)
result
##方差未知时，需要用到t检验统计量
##画一个正态分布和t分布的比较图
curve(dnorm(x),from=-5,to=5,ylim=c(0,0.45),ylab="",col="blue")
par(new=TRUE)
curve(dt(x,1),from=-5,to=5,ylim=c(0,0.45),ylab="",lty=2,col="red")
par(new=TRUE)
curve(dt(x,3),from=-5,to=5,ylim=c(0,0.45),ylab="",lty=3)
text.legend=c("dnorm","dt(1)","dt(3)")
legend("topright",legend=text.legend,lty=c(1,2,3),col=c("blue","red","black"))
##t检验
t.test(x,mu=7)
总体方差的区间估计
##编写方差区间估计的函数
chisq.var.test<-function(x,alpha){
options(digits=4)
result<-list()
n<-length(x)
v<-var(x)
result$conf.int.var<-c(
(n-1)*v/qchisq(alpha/2,df=n-1,lower.tail=F),
(n-1)*v/qchisq(alpha/2,df=n-1,lower.tail=T))
result$conf.int.se<-sqrt(result$conf.int.var)
result
}
chisq.var.test(x,0.05) ##开方后才是置信区间
2.1.3 双总体均值差的估计
根据正态分布的可加性得知，如果有
X
1
∼
N
(
μ
1
,
σ
1
2
)
X_1\sim N(\mu_1,\sigma_1^2)
X1∼N(μ1,σ12)和
X
2
∼
σ
2
2
X_2\sim \sigma_2^2
X2∼σ22，则有：
a
X
1
+
b
X
2
∼
N
(
a
μ
1
+
b
μ
2
,
a
2
σ
1
2
+
b
2
σ
2
2
)
aX_1+bX_2\sim N(a\mu_1+b\mu_2,a^2\sigma_1^2+b^2\sigma_2^2)
aX1+bX2∼N(aμ1+bμ2,a2σ12+b2σ22)，
当a=1,b=-1时，有：
X
1
−
X
2
∼
N
(
μ
1
−
μ
2
,
σ
1
2
+
σ
2
2
)
X_1-X_2\sim N(\mu_1-\mu_2,\sigma_1^2+\sigma_2^2)
X1−X2∼N(μ1−μ2,σ12+σ22)，
x
ˉ
1
−
x
ˉ
2
∼
N
(
μ
1
−
μ
2
,
σ
1
2
/
n
1
+
σ
2
2
/
n
2
)
\bar{x}_1-\bar{x}_2\sim N(\mu_1-\mu_2,{\sigma_1^2}/n_1+{\sigma_2^2}/n_2)
xˉ1−xˉ2∼N(μ1−μ2,σ12/n1+σ22/n2).
（这里和原文有出入，原文给的是开方后的标准差）。
-当方差已知时-
(
x
ˉ
1
−
x
ˉ
2
−
z
n
/
2
σ
1
2
n
1
+
σ
2
2
n
2
,
x
ˉ
1
−
x
ˉ
2
+
z
n
/
2
σ
1
2
n
1
+
σ
2
2
n
2
)
(\bar{x}_1-\bar{x}_2-z_{n/2}\sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}},\bar{x}_1-\bar{x}_2+z_{n/2}\sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}})
(xˉ1−xˉ2−zn/2n1σ12+n2σ22,xˉ1−xˉ2+zn/2n1σ12+n2σ22)
-方差未知时-
(
x
ˉ
1
−
x
ˉ
2
−
z
n
/
2
s
1
2
n
1
+
s
2
2
n
2
,
x
ˉ
1
−
x
ˉ
2
+
z
n
/
2
s
1
2
n
1
+
s
2
2
n
2
)
(\bar{x}_1-\bar{x}_2-z_{n/2}\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}},\bar{x}_1-\bar{x}_2+z_{n/2}\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}})
(xˉ1−xˉ2−zn/2n1s12+n2s22,xˉ1−xˉ2+zn/2n1s12+n2s22)
chicks<-data.frame(feed=rep(c(1,2),times=c(3,6)),weight_gain=c(42,68,85,42,97,81,95,61,103))
tapply(chicks$weight_gain,chicks$feed,mean)
tapply(chicks$weight_gain,chicks$feed,sd)
##t-test
t.test(weight_gain~feed,data=chicks,var.equal=TRUE)
##t.test有参数选择paired，默认FALSE为独立样本，TRUE为配对样本；参数conf.level默认为0.95的置信水平，var.equal默认为FALSE，TRUE表示两个总体具有相同方差。
##配对样本
Feed.1<-c(44,55,68,85,90,97)
Feed.2<-c(42,61,81,95,97,103)
t.test(Feed.2,Feed.1,paired=T)
##等价于
diff=Feed.2-Feed.1
t.test(diff)
2.1.4 双总体比例差估计
p
^
1
−
p
^
2
∼
N
(
p
1
−
p
2
,
p
1
(
1
−
p
1
)
n
1
+
p
2
(
1
−
p
2
)
n
2
)
\hat{p}_1-\hat{p}_2\sim N(p_1-p_2,\frac{p_1(1-p_1)}{n_1}+\frac{p_2(1-p_2)}{n_2})
p^1−p^2∼N(p1−p2,n1p1(1−p1)+n2p2(1−p2))
-比例未知时-
p
^
1
−
p
^
2
±
z
α
/
2
p
^
1
(
1
−
p
^
1
)
n
1
+
p
^
2
(
1
−
p
^
2
)
n
2
)
\hat{p}_1-\hat{p}_2 \pm z_{\alpha /2}\sqrt{\frac{\hat{p}_1(1-\hat{p}_1)}{n_1}+\frac{\hat{p}_2(1-\hat{p}_2)}{n_2})}
p^1−p^2±zα/2n1p^1(1−p^1)+n2p^2(1−p^2))
prop.test(x=c(225,128),n=c(500,400),correct=F)##corrrect表示是否需要连性续修正
2.2.3 均值检验
计算p值
##双尾检验
2*pt(-2.93,16,lower.tail=T)
##单位检验
pt(-2.93,16)
##等价于
t.test(pH,mu=7,alternative="less")
##双尾的临界值
qt(0.025,16);qt(0.975,16)
2.3 极大似然估计
极大似然估计的基本思想：设总体含有待估参数θ，θ可以取很多值，，所以要在很多可能取值中选择一个使样本观测值出现的概率最大的θ值，记为 θ ^ \hat{\theta} θ^，并称其为θ的极大似然估计。
2.3.2 求极大似然估计的方法
当函数关于参数可导时，可求导获得极大值对应的参数值，也为了求导方便，常会取对数后根据“费马定理”求偏导取零得到似然方程，求解。参数的极大似然估计并不能保证无偏性，极大似然估计还有一个“不变原则”，即参数的函数，参数的似然估计对应的函数也是对应函数的极大似然估计。
2.3.3 极大似然估计应用举例
先确定似然函数的表达式，然后用R中的极值求解函数。在单参数情况下，例如：
f<-function(lambda){
logL = n*log(lambda)-lambda*sum(x)
return(logL)
} ##这是一个指数分布的似然函数表达式
x = c(518,612,713,388,434)
n = length(x)
duration<-optimize(f,c(0,1),maximum=TRUE)
## interval指定区间内搜索函数f的极值，也可用lower和upper来控制，参数maximum=FALSE表示求极小值，TRUE表示极大值
duration ##得到参数值lambda的MLE
1/duration$maximum ## 由指数函数的期望为1/lambda 得到
###另一种方法和案例
Library(MASS)
Attach(geyser) ##需要的数据集
Hist(waiting,freq=FALSE,col=”wheat”)
Lines(density(waiting),col=’red’,lwd=2)
##因为绘制的直方图表示图形有两个峰，像两个分布叠在一起的，所以推断分布是两个正态分布的混合，此时需要估计的参数有5个
LL<-function(params,data){
T1<- suppressWarnings(dnorm(data,params[2],params[3]))
T2<- suppressWarnings(dnorm(data,params[4],params[5]))
ll<- sum(log(params[1]*t1+(1-params[1])*t2))
Return(ll)
}
##suppressWarnings()函数用于处理在迭代过程中可能遇到的无效值信息
##R中的maxLik程序包用于进行极大似然估计，其中的method参数可以选择不同的数值求解方法。
Library(“maxLik”)
Mle<- maxLik(logLik = LL,start = c(0.5,50,10,80,10),data=waiting)
Mle
##用图形看一下拟合效果
A<-mle$estimate[1]
Mu1<- mle$estimate[2];s1<- mle$estimate[3]
Mu2<- mle$estimate[4];s2<- mle$estimate[5]
X<-seq(40,120,length=100)
F<-a*dnorm(X,mu1,s1)+(1-a)*dnorm(X,mu2,s2)
Hist(waiting,freq = FALSE,col = “wheat”)
Lines(density(waiting),col=’red’,lty=2)
Lines(X,f,col = “blue”)
Text.legend=c(“Density Line”,”Max Likelihood”)
Legend(“topright”,legend = text.legend,lty= c(2,1),col = c(“red”,”blue”))
5599
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
